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FOREWORD 

This report presents the results of work performed by- 
Lockheed' s Huntsville Research & Engineering Center while 
under subcontract to Northrop Nortronics (NSL PO 5-09287) 
for the Aero-Astrodynamics Laboratory of Marshall Space 
Flight Center (MSFG), Contract NAS8-20082. This task was 
conducted in response to the requirement of Appendix E-1, 
Schedule Order No. E-86. 
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SUMMARY 

This report represents the development of a computer pro- 
gram for a preliminary analysis of the relative motion of a "free 
flying" experiment module in the vicinity of a Space Station under 
the perturbative effects of drag and earth oblateness. A listing 
of a computer program developed for determining the relative 
motion of a module utilizing the Cowell procedure is presented, 
as well as instructions for its use. 
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Section 1 
INTRODUCTION 

In this decade, large, semi -permanent, manned space stations will be 
latinched. These stations will provide the facilities to study and understand 
the nature of space as well as the bases for continuously observing the earth 
and its atmosphere. Experiment modules containing laboratory facilities will 
operate either attached to the stations or detached or "free flying, " depending 
on requirements of the experiments. 

This study was undertaken to develop a computer program to analyze 
the relative motion of an experiment module and a space station as they travel 
in orbit. The program considers a specific case in which the module operates 
in a "free flying" mode near the space station. Program capability is not 
limited to this specific application, however, as the program can be used 
in any situation in which the relative motion of two vehicles in nearly the same 
orbit is desired. For example, "booster-spacecraft" separations can be , 
examined. 

In developing the computer program, two approaches for examining re- 
lative motion appear; (1) a simplified approach, in which only two-body or 
Keplerian motion of the module and the station are considered, and (2) a more 
realistic approach in which are considered deviations in the motion of both 
vehicles due to the atmosphere and shape of the central body and external 
forces (such as those of the gravity of the sun and moon, and the sun's radia- 
tion pressure). The program as developed in the study contains only the per- 
turbative forces of the earth's shape and its atmosphere. A simple two-body 
(central force field) relationship could depict the motion of both vehicles, and 
integrating the force equations would lead to, in both instances, simple ellipti- 
cal orbits, the planes of which are fixed in inertial space. The real earth is 
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not spherical, however, and it does have an atmosphere, both of which cause 
perturbations to the Keplerian motion. 

A preliminary on-orbit sequence for the station and module which has 
been proposed is to: {1) detach the module from the station, (2) use a propul- 
sive maneuver to achieve a higher orbit and (3) circularize the module's orbit 
at some predetermined height above the station. The module is then in its 
"stationkeeping" position. In this mode, the module is in an orbit nearly 
identical to that of the station, differing only in height. It will be assumed in 
the analysis that the module has been placed in the "stationkeeping position. " 
Referring to Fig. 1, this position is shown as location A. Under the combined 
action of drag and oblateness, the gross motion of the module, relative to the 
station, is depicted in Fig, 1. Due to a larger semi-major axis, which results 
in a slower angular rate, the module initially falls behind the station. The 
larger area-to-mass ratio of the module results in a greater drag force on the 
module than on the station, resulting in loss of altitude by the module. The 
above series of events causes. the module to move from position A to position 
B. As the module continues to lose altitude, it reaches position C. At this 
point, which is the maximum recession distance, the module and the station 
are at the same altitude. As it continues to lose altitude, the velocity of the 
module relative to the station increases; it moves to position D, and finally 
catches up to the station (position E). The ..module will pass the station unless 
some maneuver is initiated to return it to its initial position (position A). 

In subsequent sections of the document are derived, the coordinates of 
the module relative to the station, perturbation techniques applicable to the 
program, and a detailed description of the station -module program. 
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Section 2 

COORDINATE SYSTEM FOR EXPRESSING RELATIVE MOTION 

A coordinate system which has its origin at the station and moves with the 
station is used. 


In Fig. 2, the positive Yl-axis is along the radius vector to the station from 
the center of the earth pointing away from earth, the positive Zl-axis is in the 
direction of the angular momentum vector of the station's orbit and the positive 
XI -axis is in a direction such as to form a right handed coordinate system. 


Referring to Fig. 3, assume unit vectors i, fe, having the directions of 
the positive X, Y, Z, axes of a three-dimensional rectangular earth- centered 
system. Assume the position and velocity vectors of both the station and the 
module are known in this system. Thus given. 


R = 
s 


AAA 

X i + Y j + Z^k 

S S'* s 


- A * A « A 

V =Xi+Yj+Zk 
s S S"^ s 


— A /V A 

R =Xi+Yj+Zk 
m m m** m 


m 


• /N > A * A 

= X i + Y 3 + Z k 

TYn TYT*' -m 


where R and V are the radius vector and velocity vector of the space station 
s s 

and R , V are the radius vector and velocity vector of the module, 
m m • 


The angular momentum vector of the station orbit is given by 


H 


= R X V 
s s 
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The distance 7. of the module above or below the plane of the station is given by 


where 
thus 

The' magnitude of the projection of the module's position vector onto the 
station's, plane is given by 

R' = |R cos (90 -cp) ~ |R I sin<^ 
m 1 I 

Taking R^ x results in a vector Q which lies in the station's plane and is 
perpendicular to 

The angle 0 between the Q vector and the station's position vector R 

3 

is given by 

COS0 = Rg • q/|rJ |q| 

Thus the module's x-distance and y-distance can be found by 

X = R' sin (90 - 9 ) = R' cos9 
m m 

y = R' cos (90 - 9) - R I = R' sin0 - |r 
^ m ' ' s| m Is 

As a matter of convention, all positive resiolts will indicate the module to be 
behind station, higher than or above plane of station; for example, refer to 
Fig, 3, where the x, y, z distances are shown, x is negative, y is negative, 
and z is positive. The derivations are not restricted to coplanar orbits. 
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Section 3 

METHODS OF COMPUTING RELATIVE MOTION 

In Section 2, the coordinates of the naodule relative to the station were 
derived referenced to a rectangular coordinate system. A method is now 
needed to compute these coordinates continuously, including drag and oblate- 
ness perturbations. 

Two basic classes of methods or perturbation are available: "special 
perturbations" and "general perturbations. " In "special perturbations, " ac- 
celerations of the disturbed body are integrated by using numerical techniques. 
Consequently, these methods generate a particular orbit for a particular dis- 
turbed body, for particular initial conditions. The methods are ideally suited 
for calculating orbits having limited duration. Utilizing a step-by-step pro- 
cess, the perturbed orbit is continuously determined. The methods of special 
perturbations are usually classified according to the formulation of the equa- 
tions to be integrated. Two examples of these formulations are Cowell's 
method and Encke's method. The main drawback of special perturbations is 
that errors accumulate from truncation and roundoff. Truncation error results 
from the difference between the exact solution of the difference equations 
which approximate the differential equations themselves, whereas the round- 
off error results from the difference between the computed and the exact solu- 
tions of the difference equations. In numerical integration these errors are 
difficult to control. 

General' perturbations are concerned with analytical methods in which 
the accelerations are expanded into series and integrated term by term. These 
methods result in solutions to the equations of motion in the form of symbolic 
formulas which express the sought -for quantities as explicit functions of either 
time, constants of the problem or constants of integration. Examples of these 
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formulations are the "variation of coordinates" and "variation of parameters" 
(usually orbital). The methods of general perturbations are ideally suited for 
the prediction of orbits extending over many periods. The main disadvantage 
of such methods is that most contain terms for the effects of the disturbing 
potential but do not include the effects of drag, or if drag is included, it is a 
simplified drag model. Atmospheric density is generally expressed only as 
a simple exponential function of altitude and in some formulations, is applied 
to the drag equation only at perigee. 

A "variation of parameter" (general perturbations) formulation was selected 
from Ref. 1 and applied to a representative station -module example case. The 
equations were quite similar to those of Kozai (Ref. 2). Singularities in the 
equations occur for equatorial orbits, circular orbits, and orbits at the critical 
inclination. The results from the test case indicated that for perturbed motion, 
where information from point to point along the perturbed orbits is needed (time 
intervals of five minutes were used), general perturbation methods are not 
accurate enough for studying the relative motion between two vehicles in nearly 
the same orbit. Results from the test case are discussed in Section 5. Equations 
are presented in Appendix B. 

The Cowell (special perturbation) formulation was selected to generate 
the geocentric rectangular coordinates of the station and the module. From 
these coordinates, the relative coordinates of the module can be determined 
as outlined in Section 2. 

A numerical integration scheme is usedto integrate the total acceleration equ- 
ations in the Cowell formulation. The method is straightforward, and makes no 
distinction between the disturbing accelerations and the two-body (central body) 
accelerations. As a result, many significant figures must be carried in a manner 
that the disturbing accelerations are not overshadowed by the central body accelera- 
tion in the numerical integration procedure. A small integration step size (30 
seconds for this analysis) should be used to minimize the truncation error. 
However, with a small integration step size and a large number of steps, the 
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influence of roxmd-off error will be prominent. Thus, this procedure is re- 
stricted to calculation of orbits having a duration of only a few days. 

The numerical integration of the equations of the Cowell formulation is 
performed by a fourth order Runge-Kutta numerical procedure. The equations 
of the Cowell formulation are given in Appendix A. 
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Section 4 

THE STATION-MODULE PROGRAM 

4.1 DISCUSSION 

A program incorporating the Cowell formulation to compute the relative 
motion of the station and module was developed in double precision for the 
Univac 1108 (Exec 8) computer system. The two-body relative motion, as well 
as the perturbed relative motion, is determined in the program. The effect 
of the gravitational zonal harmonics through the fourth (J^) are considered. The 
density values for the drag perturbations are computed by the MSEC Modified 
Jacchia Model Atmosphere (1967) which is recommended in Ref. 3. 

4.2 INPUT 

Initial input to the Station- Module (STA-MOD) program needed to execute 
the program successfully are the semi-major axis, eccentricity, inclination, 
ascending node, argument of perigee and true anomaly of the station and the 
module. This information is input on the first two input cards, respectively. 
The elements are then transformed to position and velocity coordinates in a 
geocentric rectangular coordinate system for use in the Cowell scheme. The 
ballistic coefficients (Cj^A/m) of the station and module respectively, needed 
for use in the drag calculations, are input on the third card. 

Control- of the integration step size in the Cowell scheme, cutoff time, 
and print time is inserted on input Card 4. 

The Modified Julian Date at which the initial orbital elements were 
determined is input on Card 5. This date is necessary for use in the Jacchia 
density model to determine the semi-annual variation of density. 
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Atmospheric density is affected to a great extent by variable solar 
activity and geomagnetic activity. Consequently it must be corrected to 
account for these phenomena. The index of solar activity, the 10.7 cm 
decimetric flux, is input as a function of the year to be utilized in the density 
calculation. The total number of solar flux values and corresponding year 
values are input on Card 6, whereas their values are input on Card 7. The 
geozhagnetic activity index, Ap, is an average value determined within the 
program based on the value of the solar flux. 

Plots of the perturbed orbital elements of the station and module, the 
perturbed coordinates of the module relative to the station, the two-body re- 
lative coordinates of the module, and the deviation of module coordinates 
from two-body behavior are available. 

Table 1 gives the input cards necessary to execute the program; Table 
2 illustrates a 1108 run request with instructions for plots. 

A complete program listing is given in Appendix C. 

4.3 OUTPUT 

Table 3 gives an example of the output from the STA-MOD program. In 
the first block of data the orbital elements of the station PNUIS, true anomaly, 
AIPS, semi-major axis, EIPS, eccentricity, FINCPS, inclination, CAPWS, 
ascending node, SMAWS, argument of perigee, MEANPS, mean anomaly are 
given. The time in minutes, TTIMEM, and time in days,TTIMED, are also 
shown. 

In the second block of data are given the corresponding elements for the 
module: PNUIM, AIPM, EIPM, FINCPM, CAPWM, SMAWM, and MEANPM. 
Time for the module and station is the same. 

In the third and final block for a given time is the relative distance, in 
kilometers, of the module from the station in the x-direction, DEL TAX; the 
y-direction, DELTAY; and the Z -direction, DELTAZ. The deviation in the 
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relative motion of the module from two -body relative motion is shown next. 
The deviation is given with respect to the three coordinate distances. They 
are DEVX, DEVY, and DEVZ. A new block of data begins after this block. 

Computer plots of the orbital elements and relative motion plots of the 
modules are obtained as an output. Examples of the plots are given in 
Section 5. 

The output shown in Section 5 (Table 3) was generated by the first two 
input data cards shown in the program listing. 
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Section 5 

RESULTS AND CONCLUSIONS 

To illustrate the functional ability of the program, results of representa- 
tive station module cases are presented. 

For the first case, the effects of drag were not considered and a non- 
circular orbit was chosen. An orbit of this type was selected so that results 
of using Koelle's general perturbation equations and the Cowell .special per- 
turbation technique could be compared. There is no provision in Koelle's 
equations for drag effects and there is a breakdown in the computations for 
circular orbits. 

The initial orbital elements are; 



station 

Module 

a 

7642.45 

7642.655 

e 

0.1 

0.1 

i 

55 deg 

55 deg 


0 deg 

0 deg 

CO 

0 deg 

0 deg 

e 

0 deg 

0 deg 


The module is initially .186 km above the station. The stations initial 
perigee altitude is 500 km (270 n,mi.). 

Figures 4 through Fig. 24 depict the results for the above case in which 
the Cowell special perturbation formulation was used. The time period used 
(900 minutes) was completely arbitrary. Figures 4, 5, and 6 show the module's 
X, y and z relative distances versus time, respectively. Figure 7 shows the 
y-relative distance versus the x-relative distance. This is the view of the 
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module in the plane of the station as viewed from the station. In Figs. 4, 5 and 
6 the short period variations are quite evident. 

Figures 8, 9 and 10 depict the variation from two-body behavior of the 
module's relative coordinates. It can be seen that, for this orbit and over the 
time period given, this variation grows to one on the order of 100 meters for 
the x~relative coordinate and y-relative coordinate, but of a magnitude of 10 
meters for the z coordinates. Since for two-body motion, no motion exists out 
of plane. Figs. 10 and 6 are identical. Figures 11 through 24 show the perturbed 
orbital elements of the station and the module. 

Figures 25 through 45 represent results for the same case as above using 
Koelle's general perturbation equations. Figures 25, 26 and 27 show the modules 
X, y, and z relative distances versus time. Figure 7 gives a view of the module's 
motion in the plane of the station. When these plots are compared with the 
corresponding plots generated by using Cowell's formulation, good agreement 
is found between the y and z relative distances. In Fig. 25 the x-relative coor- 
dinate (Koelle's equations) appears to have an additional periodic variation super- 
posed on the "known" short period variation. Many procedur es w ere instituted 
in an effort to remove this additional wiggle, but to no avail. Figure 28 and 
Fig. 7, depicting the motion in the station plane, do not agree, because Koelle's 
x-distance does not match Cowell's x-distance values. 

The deviations from two-body behavior (Figs, 29, 30 and 31) agree fairly 
well with Cowell corresponding plots (Figs. 8, 9, 10) only with the deviation from 
two-body behavior in the z-relative coordinates. The Koelle's plots (Fig. 29 
and 30) do not agree in form or magnitude. Figures 32 through 45 depict the 
perturbed orbital elements of the station and module. It should be remembered 
that only the perturbative effects of oblateness were considered in the orbit 
discussed above. 

As an additional example in which drag effects are included, the initial 
conditions of the station and module (Table 4) were considered. The station 
is 500 km (270 n,mi-) above the earth and the module is ,152 km (500 ft) above the 
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station. The ballistic parameter of the module (Cj^A/M) is 2.5 times that of 
the station. The drag coefficient has been taken to be 2.2. Since the date of 
the initial conditions is a future one, the solar activity index, FTENB, will 
be a predicted value based on mean of past values. Results to be presented 
are from Cowell formulation. 


Table 4 

EXAMPLE CASE 2 


Parameter 

Semi-major Axis, a 
Eccentricity, e 
Inclination, i 
Ascending Node, 
Argument of Perigee, o) 
True Anomaly, v 
Ballistic Parameter, C 


Station 

6878.566 km 

0.0 

30.0 deg 
0.0 deg 
0.0 deg 
0.0 deg 

A/m 0.0082 mVkg 


DATE: MAY 1, 1980 


Module 

6878.7084 km 
0.0 

30.0 deg 
0.0 deg 
0.0 deg 
0.0 deg „ 
0.0205 mVkg 


The example is a "loop case" where the module falls behind the station 
and subsequently catches up. Figure 46 depicts behavior under the example 
conditions, and shows the module leading the station after approximately 
3750 minutes. Figure 47. show the y-relative coordinate. The magnitude of 
the fluctuation about a mean value appears to increase as the module falls 
below the station. The z-relative coordinate. Fig. 48, tend to fluctuate about 
the station's plane, returning to this plane in the same terms as it takes the 
module to complete its loop. Figure 49 illustrates the loop. 

Note that this example was taken, in total, from Ref. 4. The Ref. 4 
analysis indicated that the maximum recession distance would be (27.8 km) 

(15n. mi. ) and that the time in the "far-out loop" would be 4.63 days. The 1959 
ARDC Density Module, with corrections for solar activity, was used to com- 
pute density values. The MSFC Modified Jacchia Model Atmosphere (1967) was 
used in this analysis. 
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Figure 50 shows the y-relative distance vs. x-relative distance results 
utilizing the Earth Orbital Decay program (Ref. 5). This program utilizes a 
first-order variation of parameters technique, in which, the short period 
variations have been averaged out^ The result is that which would be obtained 
if a mean line were drawn through the results in Fig. 49. The time for the 
"far-out loop" from this procedure was 2.75 days. The Jacchia model was 
used in this procedure also. 

The results of the program indicate that the program can be used to 
investigate the behavior of an experiment module or any other vehicle relative 
to another moving vehicle. The disparity between the general perturbation 
scheme and the special perturbation scheme should be resolved. In addition, 
the addition of some type of propulsive capability would be extremely helpful 
for program flexibility. 

Accuracies or inaccuracies which could be incurred when utilizing the 
program were not determined. 
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Table 1 

STA-MOD PROGRAM INPUT 


Input 

Card Program Symbol 

No. 

1 AiPS, EIPS, FINCPS, 
CAPWS, SMAWS, 
PNUIS 


2 AIPM. EIPH 

FINCPM, CAPWM, 
SMAWM, PNUIM 


3 CDAS, CDAM 


4 DT, TCUT, NP 


5 XJD 


6 K 


7 FTENB 


Definition 


Semi-major axis (km), eccentricity, 
inclination (deg), ascending node (deg), 
argument of perigee (deg), and true 
anonaaly of station (6E12. 8) 

Semi-major axis (km), eccentricity, 
inclination (deg), ascending node (deg), 
argument of pCrigee (deg), and true 
anomaly (deg) of module (6E12, 8) 

Ballistic coeffi^ent of station and 
module (meter /kg) (2E12. 8) 

Integration step size (sec) for Cowell 
method, cutoff time (hr), number of 
DT’s (integer) per print interval (print- 
out will occur every DT x NP seconds) 

(2E12. 8, 13) 

Modified Julian date at which initial 
orbital elements for the station and 
module were given (E12. 8) 

An integer which specifies the total 
number of values to be read on card 
no. 7 (13) 

Table of 81-da^^mean values of the 10. . 7 cm, 
solar flux (10 watt s/m' -/cyl/ sec). For 

future flights, predicted values are input. 
The values are loaded in the order FTENB, 
decimal year, FTENB, decimal year, etc. 
up to 100 values and the corresponding 
year may be Loaded 
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Table 2 

UNIVAC 1108 INSTRUCTIONS FOR PLOTS 


1108 RUN REQUEST 8. INSTRUCTIONS 


NAME (LAST 1. INITM 

IBtV . z 

iL) 

r. 

OP 

CODE 

/a 

JOBS PROD n 

<ffo5VC5 

BINU BLDG# 

RMS 

t 1 

RUN-ID 

3ut^^s'/z, 

RUN , 
/ OF / 

EST. CPU RUN TIME 
HRS. _^MINS.: 

CORE SIZE PUNCH 

52/^ 

iS COMPILE 
JB EXECUTE 
O SORT 

Q LANGUAGE MAX. 

S: exec VIII 

□ SPECIAL FORMS 
TYPE COPIES 


DOES THIS JOB HAVE A RESTART PROCEDURE? □ YES CJ NO 



ST>9 rVOAJ — /^oDU^ ^ 

□ OVER 


MICUO FILM 

COPIES 

COPY FLO 

OPER. miT. 

(FILES UFRAMES 

P 

F 



SEO.# 

1 

Pvl 

/ 


/ 




OPERATOR COMMENTSi d SEE TECH. C3 SEE OPeR. 


MBFC • Pom aotj (Pov Augoit 19B9) 


□ OVER 


17 

LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D1 62646 


FIGURES 

Page 

1 Relative Motion of Module as Seen From Space Station 19 

2 Relative Coordinate System 20 

3 Rectangular Earth- Centered System 21 

4 Modules x -Relative Position {km) vs Time (min) 22 

5 Modules y -Relative Position (km) vs Time (min) 23 

6 Modules z -Relative Position (km) vs Time (min) 24 

7 Modules y -Relative Position (km) vs Modules x-Relative 

Position (km) — Motion in Station's Plane 25 

8 Deviation in Modules x-Relative Position from Two-Body 

Relative x-Position (meters) vs Time (min) 26 

9 Deviation in Modules y-Relative Position from Two-Body 

Relative y-Position (meters) vs Time (min) 27 

10 Deviation in Modules z-Relative Position from Two-Body 

Relative z- Position (meters) vs Time (min) 28 

11 Stations Semi-Major Axis (km) vs Time (min) 29 

12 Stations Eccentricity vs Time (min) 30 

13 Stations Inclination (deg) vs Time (min) 31 

14 Stations Ascending Node (deg) vs Time (min) 32 

15 Stations Argument of Perigee (deg) vs Time (min) 33 

16 Stations True Anomaly (deg) vs Time (min) 34 

17 Stations Mean Anomaly (deg) vs Time (min) 35 

18 Modules Semi-Major Axis (km) vs Time (min) 36 

19 Modules Eccentricity vs Time (min) 37 

20 Modules Inclination (deg) vs Time (min) 38 

21 Modules Ascending Node (deg) vs Time (min) 39 

22 Modules Argument of Perigee (deg) vs Time (min) 40 

23 Modules True Anomaly (deg) vs Time (min) 41 

24 Modules Mean Anomaly (deg) vs Time (min) 42 

25 Modules x-Relative Position (km) vs Time (min) 

(Koelle' s Equations) 43 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D1 62646 


Figures (Continued.) 

Page 

26 Modules y-Relative Position (km) vs Time (min) 

(Koelle's Equations) 44 

27 Modules z -Relative Position (km) vs Time (min) 

(Koelle's Equations) 45 

28 Modules y-Relative Position (km) vs Modules z -Relative 

Position (km) ~ Motion in Stations Plane (Koelld s Equations) 46 

29 Deviation in Modules x-Relative Position from Two-Body 

Relative x-Position (meters) vs Time (min) (Koelle's 
Equations) 47 

30 Deviation in Modules y-Relative Position from Two-Body 
Relative y-Position (meters) vs Time (min) (Koelle's 

Equations) 48 

31 Deviation in Modules z-Relative Position from Two-Body 
Relative z-Position (meters) vs Time (min) (Koelle's 

Equations) 49 

32 Stations Semi-Major Axis (km) vs Time (min) (Koelle's 

Equations) 50 

33 Stations Eccentricity vs Time (min) (Koelle's Equations) 51 

34 Stations Inclination (deg) vs Time (min) (Koelle's Equations) 52 

35 Stations Ascending Node (deg) vs Time (min) (Koelle's 

Equations) 53 

36 Stations Argument of Perigee (deg) vs Time (min) (Koelle's 

Equations) 54 

37 Stations True Anomaly (deg) vs Time (min) (Koelle's 

Equations) 55 

38 Stations Mean Anomaly (deg) vs Time (sec) (Koelle's 

Equations) 56 

39 Modules Semi-Major Axis (km) vs Time (min) (Koelle's 

Equations) 57 

40 Modules Eccentricity vs Time (min) (Koelle's Equations) 58 

41 Modules Inclination (deg) vs Time (min) (Koelle's Equations) 59 

42 Modules Ascending Node (deg) vs Time (min) (Koelle's 

Equations) 60 

43 Modules Argument of Perigee (deg) vs Time (min) 

(Koelle's Equations) 61 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D162646 


Figures {Continued) 

Page 

44 Modules True Anomaly (deg) vs Time ,(min) (Koelle's 

Equations) 6Z 

45 Modules Mean Ajiomaly (deg) vs Time (min) (Koelle's 

Equations) 6 3 

46 Modules x- Relative Position (km) vs Time (min) — Loop Case 64 

47 Modules y-Relative Position (km) vs Time (min) — Loop Case 65 

48 Modules z-Relative Ppsition (km) vs Time (min) — Loop Case 66 

49 Modules y-Relative Position (km) vs Modules x-Relative 

Position (km) — Motion in Stations Plane 67 

50 Y- Component vs X- Component (Earth Orbital Lifetime 

Deck) 68 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



Table 3 

PROGRAM OUTPUT 
•••••DESCRIPTION OF DATA OUTPUT^^*«* 

station true anomaut idegsi • - - - — 

station SEHIoHAUOR axis (KNSI station eccentricity station- INClIHATlON (OEOSI - STATION LON« OF ASCEUOlNC NODE (DECS) 

station ARCUMENT of perigee IOEGS) STATION MEAN ANOMAUV ( PEGS ) 


station time ININS) 

station TINE lOAYS) 


• 



-HOOULC TRUE ANONaUV IOCGSI 







module SENt'NAJOR AXIS (KHS) 
NODULE ARGUMENT OF PERIGEE IOEGS) 

NODULE ECCENTRICITY 
MODULE mean ANOMALY IOEGS) 

module inclination IOEGS) 

MODULE 

LONG OF ascending NODE IDEGSI 

DELTA! OF 
DEVIATION 

nodule to station IKMS) 
OF TttO-BOOY-X 

oeltay of 
DEVIATION 

MOO TO STAT IKM) 
OF TWObBODY-Y 

DELTA! OF 
DEVIATION 

HOO To STAT (KM) 
OF TW0-B0DY*Z 



PNUIS ■ 
Alps - 
shams • 
ttimcn- 

• H92S76t-OS 

• 7GRS<«s6o«0<t 
-•l‘i92S7||oOS 

•OdOOOOOO 

EIPS - 
HEANPS* 

ttineo- 

.10000000*00 

•OOQOOOOO 

iOOOOOOOO 

FINCPS* 

•5SOn0008*02 

CAPRS B 

.OOOOOOOO 

PNUIN ■ 
AlPH • 
SMARM > 

*9S4S7«95->S* 
• 76<l26S5a^uN 
oi9SAS770B-04 

EIPH • 
KEANpm- 

1 10000000*00 
,00000000 

FINCPHb 

.550nOOOS*02 

CAPHS ■ 

•OOOOOOOO 

DELTA!" 
DEV! - 

’••3GS4zg36>|N 

-.3SS62S36-1<I 

DELTAY« 
OEVY « 

• 18519930*00 
. 18519930*00 

DELTAZ* 
DEVZ ■ 

.OOOoOOOO 

•boQooooo 



PNUIS • 

Alps ■ 

SMAWS • 
TTIMEM» 

• 19A56SA2«02 
• 74<t0<I970«0N 
•2 ZHSo2|7«66 
•SOOUQOOO+01 

eips » 
HEANPS* 

ttimed- 

.997I403I-01 

•14219194*02 

•31722222*02 

FINCPSB 

• 51994030*02 - - 

CAPRS B 

*.11386324-02 

PNUIN ■ 

*I94SS7a3*u2 







AIPH ■ 
SMAMN > 

t76<l07028»0<t 

•22NN6332«u6 

EIPH « 
HEANPM- 

•99714015-01 
• 14218539*02 

FINCPHB 

♦51994030*02 

CAPMS « 

-.11384125*02 

OELTAX. 
OEVX *• 

•96224fQ3>ol 
•942IS2I l-oi 

oeltay- 

OEVY ■ 

•18324811*00 

•18321490*00 

DELTAZ" 
DEVZ b 

•21583711*05 

*.31797803-11 



PNUIS ■ 

Alps • 

SHAMS • 
TT1MEM« 

•39tSlHA8«02 
t763S9R36^uN 
• 2t)l23ai3-»o6 
iIUO0000U*02 

EIPS ■ 
HEANPS» 
TTIMEDb 

<99152331-01 

<32527141*02 

<69111111-02 

FINCPSB 

.51984310*02 

CAPRS B 

-.81212531*02 

PNUIN • 
AlPM » 

shanh « 

• 39|N»»22+«i2 
.7636 ir97*oR 

• 20|236‘I9»u6 

EIPM « 
MEANPM" 

<99152389-01 

•32524117*02 

F INCPH* 

.51984312*02 

CAPMS >• 

-.81229381-02 

0ELTA!> 
DEVX ■ 

•18952971*00 
• IB919|86*00 

oclTay- 

DEVY b 

•17810840*00 

•17791208*00 

OELTAZb 
DEVZ B 

.3I2i<0149*D1 

-.43595404*11 


■ 

PNUIS ■ 
A I PS - 
SHANS • 
TTIMEM- 

.58138172*02 
.74311194*01 
•>13074U00>U1 
. 1SOOOOOO*02 

EIPS « 
MEANPSb 
TTIMEOb ' 

•98814209-01 

.18831589*02 

.10114447-01 

FINCPSB 

.51975792*02 

CAPRS a 

p, 

*.23009100-0) 

n 

< 

a: 


w 

M 

O 


o 
»— « 

N 

O' 

O 


PNUIH 


&6l329224>ci2 


ORIGINAL PAGS l3' 
OF POOR QUAUTB ' 



-HUNTSVILLE RESEARCH & ENGINEERING CENTER 



Fig. i — Relative Motion of Module with Respect to Space Station 
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Fig. 4 — Modiales X-Relative Position (km) vs. Time (min) 
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Fig. 5 - Modules Y-Relative Position (km) vs. Time (min) 
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Fig. 6 -Modvaes Z-Relative Position (km) vs. Time (min) 
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Fig. 8 - Deviation in Modules X-Relative Position from 

Two-Body Relative X-Position (meters) vs. Time (min) 
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Fig. 10 — Deviation in Modules Z-Relative Position from Two-Body 
Relative Z Position (meters) vs. Time (min) 
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Fig. 11 — Stations Semi-Major Axis (km) vs. Time (min) 
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Fig. 12 — Stations Eccentricity vs. Time (min) 
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Fig. 14 — Stations Ascending Node (deg) vs. 


Time (min) 
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Fig. 15 — Stations Argument of Perigee (deg) vs. Time (min) 
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Fig. 16 — Stations True Anomaly (deg) vs. Time (min) 
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Fig. 17 — Stations Mean Anomaly (deg) vs. Time (min) 
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Fig. 18 — Modules Semi-Major Axis (km) vs. Time (min) 
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Fig. 19— Modules Eccentricity vs. Time {min) 
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Fig. 21 Modules Ascending Node (deg) vs. Time (min) 
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Fig. 22 — Modules Argument of Perigee (deg) vs. Time (min) 
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Fig. 23 — Modules True Anomaly (deg) vs. Time (min) 
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Fig. 24 - Modules Mean Anomaly (deg) vs. Time (min) 


42 

LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 




LMSC-HREC D162646 



Fig, 25 — Modules X-Relative Position (km) vs. Time (min) 
(Koelle's Equations) 
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Fig. 26 — Modules Y-Relative Position (km) vs. Time (min) 
(Koelle’s Equations) 
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Fig. 27 — Modules Z-Relative Position (km) vs. Time (min) 
(Koelle's Equations) 
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Fig. 28 - Modules Y-Relative Position (km) vs. Modules X-Relative 

Position (km) - Motion in Stations Plane (Koeiles Equations) 
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Fig. 29 — Deviation in Modules X-Relative Position from Two- 
Body Relative X-Position (meters) vs. Time (min) 
(Koelle's Equations) 


47 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 


o o ^ > ‘•4 r mo I -< >-#rPio 


LMSC-HREC D162646 



Fig. 30 Deviation in Modules Y-Relative Position from Two- 
Body Relative Y Position (meters) vs. Time (min) 
(Koelle's Equations) 
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Fig. 31 - Deviation in Modules Z-Relative Position from Two-Body 
Relative Z Position (meters) vs. Time (min) 

(Koelle's Equations) 
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Fig. 32 — Stations Semi-Major Axis (km) vs. Time (min) 
(Koelle's Equations) 
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Fig. 33— Stations Eccentricity vs. Time (min) 
(Koelle' s Equations) 
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STATION IHCUHATION VS TTIH6 


Fig. 34 — Stations Inclination (deg) vs. Time (min) 
(Koelle's Equations) 
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STATION (WOE VS TIME 


Fig. 35 — Stations Ascending Node (deg) vs. Time (min) 
(Koelle's Equations) 
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station A»6 of PEhiOEE vs TIH 


Fig. 36 — Stations Argument of Perigee (deg) vs. Time (min) 
(Koelle's Equations) 
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STATION TSUE ANOMALY VS TIKE 

Fig. 37 — Stations- True Anomaly (deg) vs. Time (min) 
(Koelle's Equations) 
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STATION KSAN ASOKALf VS TIKE 


Fig. 38 — Stations Mean Anomaly (deg) vs. Time (sec) 
(Koelle's Equations) 
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XOOUUE AXIS VS TIME 

Fig. 39 — Modules Semi -Major Axis (km) vs. Time (min) 
(Koelle's Equations) 
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HCOUUe eccentricity vs TIhE 


Fig 40 - Modules Eccentricity vs. Time (mxn) 
(Koelle' s Equations) 
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HCOULE tKCLtttoTlON VS T T»t« 


Fig. 41 — Modules Inclination (deg) vs. Time (min) 
(Koelle's Equations) 
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MODULE ARt OF PERICEE VS TIME 


Fig. 43 — Modules Argument of Perigee (deg) vs. Time (min) 
(Koelle's Equations) 
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wooule true anomaly va time 


Fig. 44 — Modules True Anomaly (deg) vs. Time (min) 
(Koelle' s Equations) 
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hooule HEAN anomaly v& time 


Fig. 45 — Modules Mean Anomaly (deg) vs. Time (min) 
(Koelle's Equations) 
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Fig. 46 — Modules X-Relative Position (km.) vs. Time (min)-Loop Case 
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Fig. 47 — Modules Y-Relative Position (km) vs. Time (min;rIL.oop Case 

65 

LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 


LMSC-HREC D162646 


oo: 000 



Fig. 48 — Modules Z-Relative Position (km) vs. Time (min)-loop Case 
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OetTAX tH KH3. 


Fig. 49 — Modules Y-Relative Position' (km) vs. Modules X-Relative 
Position (km) — Motion in Stations Plane 
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Fig. 50. - Y-Component vs X-Component (Earth Orbital Lifetime Deck) 
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COWELL METHOD 
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Appendix A 


The equations of motion to be employed are 


^ - ^CBE ^2HE ^3HE ^4HE ^DRAG 


Y = Y +Y +Y +Y +Y 

'‘CBE ^ ^2HE 3HE ^4HE ^ ^DRAG 


^ = ^CBE ^2HE ^3HE ^4HE + ^DRAG 


The central body (earth) terms are 


X 


CBE 


jiX 

r3 


CBE 


= . ^ 


'CBE 


R^ 


The terms due to the second harmonic of the earth's potential are 
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2HE 


- /JXA, 
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1 - 
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'2HE 
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5Z 
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The terms due to the third harmonic of the earth's potential are 


HHXZ 


'3HE 


3HE 


'3HE 


^^4 L T2^\ 

\ ■ R^ / 


/iHYZ 




R 


-3/iHAg^^ 10Z‘ , 

■i “ • 
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(■ 
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35Z' 

3R^ 
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The terms due to the fourth harmonic of the earth’ s potential are 


X 


4HE 


A . , ^ , z"; 

R^ V R^ r"^ 


-aDY /3 


•4HE 


R 


Z^ Z^' 

— - 6 — + 9 — 
7 ° 2 ^ ^ 4 

R R / 


"4HE 


/l5 10 j . 

r7 \7 - ^ ' 1,4 


The drag perturbation terms are 


^DRAG 

^DRAG 

^DRAG 


2m 


fcl 

2m 


2m 


PV^ (X + 0>Y) 10^ ’ 
PV^ (Y - ioX) 10^ 

PV Z 10^ 
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where 

w 

•^E 

J, H, D 

where 

J 

H 

D 

P 

V 

e 

A 

m 


rotational velocity of earth’s atmosphere (rad/sec) 

mean radius of earth ellipsoid (km) 

3 2 

earth gravitational constant (km /sec ) 

second, third, fourth constants of the earth potential 
function 

3/2 ^2 = ^ 

-5/2 J 3 = +0.575 X 10"^ 

-15/4 = + 0.795 X 10"^ 

3 

atmospheric density (kg/m ) 

inertial velocity (km/sec) 

dimensionless drag coefficient 

2 

area (frontal) of vehicle (m ) 
mass of vehicle (kg) 


The equations of motion (X, Y, Z) are integrated numerically using fourth- 

order Runge-Kutta integration to establish the geocentric space-fixed 

velocities (X, Y, Z) and position (X, Y, Z). The units on all acceleration 

2 

terms are km/sec . 
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Appendix B 

KOELLE'S VARIATION OF PARAMETERS 
GENERAL PERTURBATION TECHNIQUE 
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Appendix B 

Given initial values of semi-major axis a^, eccentricity e^ , inclination 
i^, argument of perigee ascending node and mean anomaly then 

at any given time, these orbital elements are given by 


a = 


a + 
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3 2 

•j kJ A ^ 

2 2 e 


a(l - e ) 
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(continued) 
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(continued) 
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(continued) 
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35 h K- 


4-5 


It" 

sm 1 ( 


2 3-1 2 

sini - 3 e(l + e ) sin ij (8 - 9 sin i) 


u. o 3 2. 

+ 2 e cos 1 


i sini ( 


2 

ro n • 2 .. , 4 sin i 

(8 - 9 . sin i) + 


^ c ■ 2. 

4-5 sin 1 


’cosStJ 


f ^2 nt |l - tX (2 + 3 e^) ( 4-7 sin^i) 


16 V P 


+ + e^) - 5 sin^i {8 - e^)] | 




j 6 (v - M + e sinv) - 3 e sin(v + 2u)) - 3 sin( 2 v + 2co) 


e \ e cosi 


. 2 ,, 


- e sin( 3 v + 2 <J)J - '^2 1 ^ / 2(7 - 15 sin“i) 


4-5 sin i 


, c . 2 . 14-15 sin^i , ^ “^4 
+ 5 sin 1 s — + 5 Y“ 

4-5 sin i '^Z 


„ . 2 .. , c . 2 . 6 - 7 sin^i 
2(3 - 7 sin i) + 5 sin i ^ 

4 -- 5 sin i 


sin 2 o) 


1 ^3 ’‘e ,, , 5 ^ 5 /*e\ ecoti (l 0(4 + 3 e^) sin^i 

- 2 IT ^ ® ‘=°‘* + 32 jTl^) — — t 't 't: 

2 2 \ / 4 - 5 sin it 4-5 sin i 


X, 


1 ^ 8-7 sin^i (4 - 3 sin^i)J 


i (4 - 3 sin i)| cosCJ + (4 + 3 e ) 


j ^8 - 21 sin^i( 4-5 sin^i)J cosw + ~ 


2 . 2 . 
e sin 1 


■ 2 
0/0 le-• 2 .^ , - 2 . 8-9 sin i 

3(8 - 15 sin i) + 10 sin, i 2 “ 

4-5 sin i J 


cos 3 tt» I 
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M + nt 
o 




+ 8^1 - sin^i ^ 


130 25 2 ^ .„\i, 2 

_ -f48\/l-e 


«■>(-) [ 

I 


10 + 5 e' 


+ A j 2 /!;e 

^ 64 ^2\ P 


\li ^ /o 2 , . 2 . 

y 1 - e {2 - e ) sin i 


jAt 2 

128 ^41 P * ® 


F 


(8 - 40 sin\ + 35 sin^i) 


-l:h(¥^ Xfl-e" 


^[7~?'|^1 - I sin^i^ ji^i . =i 


, 1 -31 . 2 . 

+ e sm V I - sui i 


]- 


sinv + Y sin2v 


:j(~+ f ej sin(v + 2cj) - e sin(v - Zo) 


7/1 e 


3 . 


" IZ"! e” " ^ |sin(3v + 2cj) - g sin{4v + 2(o) “ Tr e sin(5v + 2oj) 


+ li <1 


3/2 . 2. 

' sin 1 


4-5 sin^i 


(14 - 15 sin^i) 


A 2 

+ 5 — ^ (6-7 sin i) 

4 


sinZ(t> + 


1 (1 -e^) 

2 jJ P J 


3/2 


sini coscj 


c J. /a „ 2,^/^ 

+ AA/_e\ (1 - e ) 


32 J 2 \ P 


(4 


^ ^ /(4-}-9e^) 8-7 sii 

- 5 sin i) ( _ ■ 


(4 + 9e^) 1 8 - 7 sin^i(4 - 3 sin^il cos u> 


J 


7 2 . 2. ,0 Q . 2., , I 

^ e sin 1 (8 - 9 sin 1 ) cos36>. 
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where 


^ 3 ’ ■^ 4 ’ '^5 



P = 


n = 
t 


earth's equatorial radius (km) 

earth's second, third, fourth and fifth geepotential 
coefficients 

+ 10,82.28 X 10"^ ' 

-2.3 X 10"^ 

-2.12x10“^ 

-0.2 X 10"^ 

true anomaly (deg) 
a(l - e^) 

mean motion of vehicle 
time 
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Appendix C 

PROGRAM LISTING 
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*RUN«//T BUTLER«4a2330«BUTLER8IN202»3«300/000 . BUTLER/LOCKHEED 

• FOR,T‘^ MAIN, main 

IMPLICIT REALMS (A-H,0-2) 

C THIS PROGRAiM DETERMINES THE RELATIVE POSITION OF AN ASTRONOMY MODULE 
C WITH respect to A SPACE STATION BY UTILIZING THE COWELL INTEGRATION 
C ROUTINE* THE PROGRAM IS CONSTRUCTED TO HANDLE PERTURBATIONS DUE TO 
C the oblate EARTH AND ATMOSPHERIC DRAG. 

COMMON/XK/NP « NDTP » DT 

O I mens ion as ( 1 4 ) , am ( 1 4 1 ♦ S a ( 1 4 ) , XM a ( 1 4 I 
■ dimension TRSA Cl 4 ) ,TBMA( 14) 

dimension SI.GMAC 14) ,ETA( 14) 

DIMENSION PLOTSt 1000)-,PLOT6I 1000) ,PL0T7( 1000) 

double PRECISION MEAnPS « MEANPM , ME'AnOS »MEANOM 

dimension PLOTC 1000) .PLOT! ( lOOO) .PLOTEC 1000) »PLOT3( 1000) 

DIMENSION PLOTXC I'^OO) ,PLOTYC 1000 ) tPLOTZC lOOC ) 

COMMON/D AT AB/XLABC 12) ,YLABl C 12) , YLAB2( 12 ) , YLAB3C 12) ♦ YLAB4 < 1 2 ) , 

1 YLABBC 12 ) , YLAB6C 1 2) , YLAB7( 12> » YLABSI 12) , YLAB9 C 12) , YLABIO C 12 ) < 
2YLAP1 1(12) »YLAB12( 12) .YLABI3C 12) iYLAB14( 12) ,YLAB1B( 12) .YLAB16C 12) « 
3YLAB17C 12) .YLABIRC 12) .YLABlPC 12) *YLAB20( 12) ,YLAB21 (12) 

Pf=-AL 

S PLOT.PLOTl ,PL0T2« PLOTS ♦PLOTX* PLOT Y.PLOTZ.XLABiYLABl , 

<SYLAB2« YLAB3* YLAB4 

real PLOTS* plots, PL0T7, YLABS « YL AB6 , YLAB7 

dimension PLOTSt lOOO) ,PL0T9( 1000) ,PL0T10C 1000) ,PL0T1 l ( lOOO) , 
1PLOT12C 1000) ,PL0 t 13( IOOO) ,PLOTl4( 1000) ,PL0T15( 1000) ,PLOT16( 1000) , 
2PLOT17C 1000) ,PL0T18( 1000) ,PL0T19( 1000) ,PL0T20( 1000) ,PL0T21 ( 1000 > 
real plots , PLOT9 • PLOT 1 O , PLOT 1 1 , PLOT I 2 * PLOT 1 3 , PLOT 1 4 , PLOT I 5 . PLOT 1 6 * 
1PLOT17,PLOT18,PLOT19,PLOT20,PLOT21 , YL ABB , YLAB9 , YL AB 1 0 , YLAB 1 1 * 
2YLAB12, YLAB13, YLAB14,YLAB15, YLAB16, YLAB17,YLAB18,yLAB19, YLAB20, 
3YLAB21 
call TDENTC9) 

0 12=1 *=^707Q<=; 

PC30 = ri,0174S329 

C READ initial ORBITAL ELEMENTS FOR STATION AND MODULE 

read ( 5,20) A I PS , E I PS , E I NCPS , CAPwS , SMAWS • PNU 1 S , A I PM , E I PM , F I NCPM , 
ICAPufM ♦ SMAWM,PNUIM 
20 FORMAT (6D1 2.0) 

C PRINT IDENTIFICATION OF OUTPUT ORBITAL PARAMETERS OF STATION AND 
C MODULE AND THE COORDINATES OF THE MODULE WITH RESPECT TO THE STATION 
WRITF< B,SO ) 
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50 F0RMAT(1H1« 40X* aeH^^-S'^'M-DESCR I PT! ON OF DATA OUTPUT*****/ 

$ P8HOFTATION TRUE ANOMALY ( DEGS ) / 

S 30H STATION SEMI-MAJOR AXIS (KMS), 7X» 20HSTATI0N ECCENTR I C I TY < 
9X» 65HSTATION INCLINATION (OEGS) STATION LONG OF ASCENDING NODE 
S (DECS) /1X» 62HSTATI0N ARGUMENT OF PERIGEE ( DECS I STATION MEAN A 
SNOMALYCDEGS) /1X» 19HSTATI0N TIME (MINS>» 1 6X » 19HSTATION TIME (DA 
YS I ) 

WRITF ( 6« 1 1 ) 

11 FORMAT! IHO, 26HM0DULE TRUE ANOMALY (DECS) / 

S 29H MODULE SEMI-MAJOR AxiS (KMS)» 8X» 19HM0DULE ECCENTR I C I TY ♦ 1 OX » 
S 68H MODULE INCLINATION (DEGS) MODULE LONG OF ASCENDING NODE ( DE 
SGS) /1X» 62HM0DULE ARGUMENT OF PERIGEE (DEGS) MODULE MEAN ANOMAL 
OiV (DFGS) • ) 

WR I TE ( ft » 1 2 ) 

IP FORMAT(1HO, 91HDELTAX OF MODULE TO STATION (KMS) DELTAY OF MOD T 
SO STAT (KM) DELTAZ OF MOD TO STAT (KM) /IXt 23HDEVIATI0N OF TWO- 
SB0DY-X» 13X» 23HDEVIATION OF TWO-BODY~Y» 6Xi 23HDEVIATION OF TWO-B 
SOOY-Z ) 

read ballistic COEFFICIENT OF STATION AND MODULE 
RFAD { 778 ) CDAS« CDAM 
77P FORM AT (2D 12. 8) 

read INTEGRATION STEP S I ZE ( SEC )♦ CUTOFF T I ME ( HR ) , I NTEGRAT I ON STEPS 
pprp ppiNT INTERVAL 
PFAD ( F . 77 ) DT « TCUT . NP 
37 FORMAT(2012.B, 13 ) 

ndtp=np 

TSTA=0.0 

tmod=o,o 

T2BFTA=0o0 
T2BMOD=C«0 

C TRAnpspoRM orbital elements to rectangular earth centered coord 

CALL TRANFM (AIPS.EIPS.fi NCPS . CAPWS . PNU I S « SMAVIS « XS . YS ♦ ZS , XDS . YDS . ZD 


•ES ) 



AS 

( 1 ) 

= XS 

AS 

(4 ) 

= YS 

AS 

(7) 

=zs 

AS 

(2) 

= XDS 
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AS<'=;)=Yr)S 

AS(8)=70S 

CALL TRANFM ( a I pm « E I pm ♦ F I NCPM « C APWM , PNU I M « SMAWM , XM « YM , ZM , XDM ♦ YDM * 2D 

AMM ) 

AM( <i > =r YM 
AM( 7) =2M 
AM( ;P ) = XOM 
AM(F^=YDM 

AM( 8)=ZDM 
8IGMA( 1 )=XS 
c;iGMA(4 )=:Y8 
F IGMA (7) =Z.S 
c;IGMA{P)=XDF 
FTGMA(=:)=YDS 
F!GMA(P)=ZOS 
FTA( 1 )=XM 
c-TA(4)=YM 
c-TA ( 7) = 71^ 
p-TA(?>=XDM 
FtA («=?) = YDM 
p=’TA CP ) = 20M 
AOF=AIPS 
FOF=FIPS 
XI0F=FINCPS 
\./OF = rAPWS 
FWOS=PMAWS 
TRU-SssONUIS 
A0M = A I <=>M 
FGM'^F I PM 
W0M=CAPWM 
XIOM=FIMCPM 
GWOM=c;MAWM 

TRijM'sPNU I M 
ttime=o. 

IPT = 0 

1 0 1 CALL COV-JELL < AS CDAS * F INCPS * CAPWS » SMAWS « MEANPS » A IPS i E I PS , PNU I S * RS * 
■»PA* TPTA ) 

IFCIPT.NE.O) GO T-^ Po°l 
NP=NDTP 
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TO OQQ? 

1 N'==n 

OOQp CONTIMUP' 

P’INCPc; = FIMCPS/RPO 
CAPW? = CAPWS/RPD 
SMAWS = =;MA<(fS/RPD 
XS=FA( 1 ) 

Y<==c;A (A ) 

7c;=PA{7) 

XDF=FA(P) 

Y0C;=:C^A(R) 

70c = c:A( fl) 

TT I TT I MF/fiO . 

TTIMFr?=TTIMF/f364 00. 

CALL COWELL( AM«CDAM«FINCPM,CAPWM,SMAWM<MEANPM » A I PM , E IPM * PNUI M » RM 
*XMA , TMOO ) 

FINCPM = FINCPM/PPD 
CAPv,(M = CAPVW/RPD 
SMAwM = SMAWM/RPO 

XM=XMA ( 1 ) 

YM=XMA(4) 

ZM=XMAC7) 

COMPUTE RELATIVE POSITIONS DELT AX« DELTAV » DELTA2 « OF MODULE TO 
station under PFDTURiSAT T VF EFFECTS 
DM=QOD -p { xm*xm+ym*ym+zm->zm ) 

HX=YS»-7DS-7S*Y0S 
HY=7S*XOS-XS*ZOS 
H7=XS*YnS-YS*XDS 
OX=YMiS-H7-7M*HY 
OY= 7M'K-HX-XM*H7 
OZ= XM*H Y- YM*HX 
HDRM = HX*XM+HY*YM+HZ-?fZM 
APH = S0RT(HX*HX+HY«-HY+HZ*HZ) 

ARQ = S0RT C OX«^OX+OY-!fQY+OZ*OZ ) 

COPH I rrHDRMZ ( A>=!H*RM ) 

SIPHI=40RT( 1 .-COPHI*COPHI ) 
xmagp=rm*c,ipht 

THETA=ACOS( (QX*XS+QY-S-YS+QZ*ZS)/( ABO*RS) ) 
nFLTAX=XMAGP-»-COS (THETA ) 

DELTAX=-DELTAX 
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DELT A Y= XMAGP-Jf-S I N ( THET A ) -RS 
DFLTA7=HDRM/ARH 

OELT,AR=SQRT(DELTAX*DELTAx+DELTAY*DELTAY+DELTA2-if-DE'LTAZ5 

C 

TF( IPT.NF.O) go to 333 
NP=N0T'= 

GO, TO 343 
•SR"? MD=n 

343 continue 

CALL TOBODY ( S I GMA - TBSA , RTBS » AOS , EOS » X I OS » WOS . SWOS . TRUS . ME AnOS . 
ITFPSTA) 

XlOS=XIOS/RPD 
SWOS=SU'OS/RPD 
WOS=WOS/PPD 
IF{IPT*NE.OI GO TO 633 

np=notp 

GO TO 643 
f.33 NP=0 
6-13 CONTINUE 

CALL TOBODY ( ETA » T8MA ♦ RTBM » AOM * EOM . X I OM , WOM » SVJOM , TRUM « ME ANOM « T2BM0D 

G WOM = s WOM /PPO 
■WOM^WOMXRPD 
XIOM=X!OM/RPD 
NP=0 

T F ( ! PT- 1 ono ) 909 « 998 , 998 

■QQQ IPT - JPT + 1 

<^9R PLOT (JPT)=TT I ME/60. 

XTBc: = TP8A ( 1 ) 

YTBS=TPSA(4I 

ZTRS=trSAC7I 

XDTRS=TRSA(?) 

YnTP8=TP8A ( R) 

7DTPS=TPSA <8) 

XTBM=TPMA<1) 

YTPM = TRMA ( <1 ) 

7TRM=tRMA (7) 

C COMPUTE relative POSITIONS DXTB .DYT8 *OZTB. OF MODULE TO STATION 

C _ FOR TWO-BODY MOTION 

H 1 = YTPS*7DTPS-ZTPS*YDTRS 
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H?=7TRS*yDTR?-XTP5^*7DTRS 

h3=xtp=^-k-ydtbs-ytrs*xdtbs 

0 1 =YTP,M*H3-7.TRM*H2 

0,P=ZTP'«1*H1 -XTBM*H3 

03=XTP!VI*HP-YTRM-S-H 1 

HRM = Hl*XTBrw'+HP*YTBM+H3*7TBM .. 

APHl = cORt<H1*H1+H2*H2+H3«H3^ ■ 

■ ABO 1=PQPT(0'1*01 +07*02+03*03 ) 

•COCH T =HRM/ { ABH 1 *RTBM ) 

CH 1= PORT (1 .-COCHI*COCHI ) 

XPTR=RTBIW*CHI 

THTB=ACOS( (Ql*XTBS+Q2*YTBS+Q3*2TBS)/(ABOl*RTES) ) 

■ nxTP=ypTR*cos cthtB) 

OXTB=-DXTB 

DYTB=XPTP*RTN<THTB) “RTBS 
07TR=HRM/ARHI 

C COMPUTE DEVIATIONS FROM TWO-BODY RELATIVE MOTION DEVX t DEVY * DEVZ 

dfvx=ofltax-dxtb 
dfvy=oeltay-dytr 
DC^V 7=nPLTA2-DZTR 
PLOTl (IPT) = DELTAX 
PLOT2 (IPT) = DELTAY 
PLOT3 (IPT) = DFLTA2 
DLOT7 ( I PT ) =DFLTAP 

C NORMAL I 2E DEVX i DEVY « DEV2 BY 10**3 THEY ARE NOW IN METERS 

PLOT^: ( I PT > =ni7VX *l.D+3 
PLOTR( IPT)=DEVY*1 .D+3’ 

PLOtT (IPT) =DEVZ* t • D+3 
RLOTR( IPT) =A I PS 
PLOTP( IPT)=FIPS 
PLOTim IPT)=:FTNCPS 
PLOTl 1 ( IPT)=CAPWS 
PLOT 1 IPT) =«.MAVJS 
PLOT 1 3 ( I PT ) =PNU I S 
PLOTlA( IPT)-MEANPS 
PLOTIFI IPT)=ATPM 
PL0T16( IPT)=EIPM 
PLOTITI IPT)=FINCPM 
PLOTIFC IPT)=CAPWM 
PLOT IP ( IPT)=SMAWM 
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PL OTpr^i IPT)=PNUfM 
PLOTPI C IPT)=MFANPf^ 

C PRINT ORBITAL ELEMENTS « T I ME » AND RELATIVE MOTION COORDINATES 

WRITE(6»60) PNUIS. A IPS* E I PS , FINCPS» CAPWS* SMAWS* MEANPS« 

% TTIMFM, TTIMFO 

60 FORMATdHO, 8HPNUIS = ♦E15.8 / 

S lx» 8HAIPS = iE15.a» 13X« 8HEIPS = ,E15.a» 6x» 8HFINCoS= ♦ 

S El5»a» 5X» 8HCAPWS = »E15.8/1X« 8HSMAWS = ,E15.8» 13X« BHMEANPS= 

S E15.8 / 9H TTIME'"= *E15.8« 13X» ShTTIMED= tElS.S) 

C 

WRITE(6»70) PNUIM# AIPM. E I PM , FINCPM^ CAPWM» SMAWM, MEANPM 
70 F0RMAT(1H0., 8HPNUIM = » E 1 5 . 8/ 

S lX« 8HAIPM = «E15,8« 13X« 8HEIPM = ,E15.8, 6X, 8HFINCPM= , 

£ EI5.8« 5X« 8HCAPMS = »E15«8 /1X« 8HSMAWM = ,E15.8» 1 3X » 

«; 6HMFAMPM= ^E1=5»F) 

C 

WRITE(6«781) OELTAX* DELTAY, DELTAZ 
781 FORMAT(1HO« 8H0ELTAX= »E15.8» 1 3X * 8HDELTAY= ,E1S.S» 6X ♦ 8HDELTAZ= 

£ ,E1E.8) 

C 

W»RITF (6« 1 OOO) DXTP» DYT8,DZTP 

1000 FORMAT(lX* 8HDEVX = «E15.8» 13Xt 8HDEVY = »E15.8» 6X i 8HDEVZ = 

£ FlF.n///) 

XNnTo=MDT° 

TT I ME=TT I ME+XNOTP-S-DT 
XY=TTTME/36O0* 

IF(XY-TCUT) 101 « 101 *102 
10? CONTINUF 

C PLOT PERTURBED ORBITAL ELEMENTS AND .RELATIVE POSITION COORDINATES 

WPITF(6*?0OP) 

POOO FORMAT! 1 HO/1 HO, 1 8HENTER I NG PLOT AREA) 

call QUIK3V f“l* IH*, XLAB, VLABl , -IPT, PLOT, PLOTl ) 

CALL QUIK3V (-1, 1 H* , XLAB, YLAB2 , -IPT, PLOT, PL0T2 > 

CALL QUIK3V (-1, IH*, XLAB, YLAB3 , -IPT, PLOT, PLOT3 > 

CALL QUIK3V(-1, 1H*» XLAB, YLAB4 , -IPT, PLOT, PLOTZ ) 
k-=o 

,VALUF = PL0T1 { 1 ) 

DO 1 «p=1,IDT 

IF(PL0T1 (KP) -VALUE) 601,600,600 

FOO 
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PLOTX ( K ) =PLOT 1 ( KP 5 
PLOTY { K ) =PLOT? ( KO ) 
1 CONTINUE 


CALL 

QUIK3V( 

-1 ♦ 


1 H*» 

YLABl 

» YLAB2* 

- K 4 

PLOTX 

4 PLOTY) 

CALL 

0UIK3V( 

-1 « 


lH*i 

XLABt 

YLA85»- 

IPT4 PLOT 4 PLOTS) 

call 

QUIK3V( 

-1 « 


iH^i-, 

XLAB. 

YLA86»- 

IPT 4 PLOT 4 PLOTS) 

call 

OU I K3V < 

-1 * 


IH*, 

XLAB« 

YLAB7«- 

IPT. plot 4 PL0T7) 

call 

OUIK3V 

<-l 


IH* 

, XLAB 

, YLAB8« 

-IPT 4 

PLOT 4 

PLOTS) 

call 

QU I K3V 

(-1 


IH* 

♦ XLAB 

« YLAB9« 

-IPT 4 

PLOT 4 

PL0T9) 

call 

QUIK3V 

(-1 


IH* 

, XLAB 

, YLABIO 

4- IPT 4 

PLOT 4 

PLOTIO) 

CALL 

OU I K3V 

<-l 


IH* 

» xlab 

♦ YLABl 1 

4 -IPT 4 

PLOT 4 

PLOT 1 1 ) 

CALL 

OU I K3V 

(-1 


IH* 

, XLAB 

« YLAB12 

4 -IPT 4 

PLOT 4 

PLOT 12) 

CALL 

QU I K3V 

(-1 


IH* 

» xlab 

♦ YLAB13 

4 -IPT 4 

PLOT 4 

PLOT 13) 

CALL 

QUIK3V 

(“1 


IH* 

« XLAB 

♦ YLAB14 

4 -IPT 4 

PLOT 4 

PL0T14 ) 

CALL 

QU I K3V 

("1 


IH* 

, XLAB 

, YLABl 5 

4 -IPT 4 

PLOT 4 

PLOT 15) 

CALL 

QU I K3V 

<-l 


IH* 

« XLAB 

4 YLAB16 

4 -IPT 4 

PLOT 4 

PLOT 16) 

CALL 

OU I K3V 

(-1 


IH* 

t XLAB 

t YLABl 7 

4-IPT4 

PLOT 4 

PLOT17) 

CALL 

QUIK3V 

(“1 


IH* 

» XLAB 

4 YLAB18 

4-IPT4 

PLOT 4 

PLOT 18) 

call 

QuIK3V 

<-l 


IH* 

< XLAS 

4 YLABIO 

4-IPT4 

PLOT 4 

PLOT19) 

call 

QU I K3V 

(-1 


IH* 

, XLAB 

4 YLAB20 

4-IPT, 

PLOT 4 

PL0T20) 

CALL QUIKSV (-1 
WRITF( pool ) 


tH* 

♦ XLAB 

4 yUAB21 

4 -IPT 4 

PLOT 4 

PLOT21 ) 


pool FOQMATMHO, 20HPL0TTING IS FINISHED) 

CALI- FNOJOB 

STOP 

END 

* FOR « I S COWELL » COwELL 

SUBROUT I NE COWELL < A » CDA « F I NCP « C APWP » SM AWP » MEANP , A I P «E I P . PNu I « 0 » Ax » 
*T ) 

IMPLICIT REALMS {A-H,0-2) 

COMMON/XK/NP ♦ NOTP » H 
dimension ■AUa)»AX{l4)»B(lF) 
nOuBLF PRECISION MpANP,YFAN 

C THIS ROUTINE INTEGRATES THE PERTURBED EQUATIONS OF MOTION USING 

C A FOURTH ORDER RUMGE-KUTTA SCHEME 

M = S 
1 M = 0 

C call MODSTA to GET INITIAL VALUES OF ACCELERATION 

CALL MODSTA ( 1 , T i A » AX I S » ECCEN ♦ ANOM . CAPW » SMAW , MEAN « F I NC , R » CDA ) 

C CALL MODSTA TO COMPUTE ORBITAL ELEMENTS FOR INITIAL AND SU3SE- 
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C QUENT TIME PERIODS ' 

CALL MODST A ( 2 » T < A » AX I S » ECCEN * ANOM i CAPW ♦ SMA W t MEAN ♦ F I NC ♦ R » COA ) 

C store POSITION* velocity. acceleration and ORBITAL ELEMENTS TO BE 

c returned to main routine 

0 =R 

AIP = AXIS 
P-IP = FCCFN 
PMUI = ANOM 
CAPiii(P=CAPW 

smawp=smaw 

FINCP = FINC 

mfanp = mean 

AX( 3 ) =A ( 1 3 
AX( 4 )=A (45 
AXf 7) =A(7) 

AX(?)=A(?) 

AX<E)=AC=:) 

AX(R)=A(8) 

C COMPUTE FIRST HALF-STEP THE FIRST TIME 

C=H/ 2 . 
kr=3*N 

n = C**2 / P»0 
F = D * 4«0 
F = C X 3.0 
T=T+C 
L = 1 

DO 2 T= 1 .X.S 
B(L)=A( I ) 

P (L+1 )=A ( I+l ) 

P(L+?)=AC I+P) 

A t I )=A( I )+C*A( l + i )+D*A( 1+2) 

A(T + n = A(I + l )+C*A( 1 + 2 ) 

2 L=L+“=^ 

C COMPUTE FIRST HALF-STEP SECOND TIME 

CALL MODST A (l.T.A.AXlS. ECCEN . ANOM , CAPVn.r . SMAw , MEAN , F I NC , R . CDA > 

L = 1 

DO 3 1 = 1 ,K .3 
B(L+3)=A( 1+2) 

A ( I+l )=B(L+ 1 )+C*A( 1 + 2 ) 

L=L+«^ 


3 



o o 


LMSC-HREC D162646 


C COMPUTE SFCOND half-stfp thf first time 

CALL MODSTAI 1 » T » A » AX I S , ECCEN i ANOM » CAPW . SMA W i MEAN » F I NC ♦ R » CDA ) 

L=1 
T = T+C 

DO 4 1=1 .K, 3 

P(L+4)=A( 1+3) 

A( r + l )=H+^A( t+3)+B{L+l ) 

A{ I)=P(L)+H^«-6{L+1 )+E*A( 1+3) 

A L=L+B 

C COMPUTE final CONDITI'ONS FOR THE TIME STFP 

CALL MODSTAI 1 » T » A » AX I S ♦ ECCEN » ANOM ♦ CAPW ,SMAW , MEAN . F I NC , R » CDA ) 

L = 1 

no T = 1 »K«3 

A ( I ) =P ( L ) +H* (B(L + 1)+F-K-(R( L+2 ) +R ( L+3 ) +B ( L+4 ) ) ) 

At I + l )=E(L+1 )+F*{B<L+2)+A( 1+2 ) +2 .•»• ( 6 ( L+3 ) +B ( L+4) ) ) 

F L^L+F 

IF(MP-NDTP) 1 »37,37 
37 return 

'FOR. IS MonSTA.MOOSTA 

SUBROUT I NE MODST a C UU ♦ T . a . Ax I S ♦ ECCEN » ANOM ♦ C APW . SMA W . MEAN . F I NC ♦ R . CD A 
*) 

IMPLICIT REAL*8 (A-H.0-7) 

COMMnN/XKXNP ♦ NOTP . DT 

THIS ROUTINE DETERMINES THE TOTAL ACCELERATION COMPONENTS AND 
ALSO COMPUTES ORbITAL ELEMENTS FROM POSITION AND VELOCITY 
LOCICAL SETI 
DOUBLE PPFCISION MEAN 

n I mens ion a ( I 4 ) 

SET I = «FALSF. 

PIE=3. 14 lEPl 
r>PR= F7 , 

A3=677F. le*^ 

W=63«6»7F4 
IF(JJ-l) 7n2.702.O0F 
702 R= SOPTt At 1 )**2+A(4)**2+Af7)**2) 

COMPUTE altitude OF VEHICLE ABOVE EARTH. S SURFACE 
AK- f 2 ) / ( A2**2 ) 

P[<=:AK'*-)f-2 

BF RHA= sORTf a f 1 ) -K-^a+A t 4 )**2 ) 
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71 

66 






0FL= ATAN( A(7)/RHA) 

PSI= ATANCAK'-K- SIN(DEL1/ cos (DEL)) 

SBI=( SIN(PSI)/ C0?(PST))**2 

SC I = SORT C C 1 ,+RK*SRI ) / ( 1 .+AK-5fSBI ) ) 

OI=A?*RCI 

ALt=P-OI 

IF( ALT-inOO, 5 65»71«71 
oHO=n , 


original 

OF POOR 


quality 


GO TO 66 
CONTIMUF 

call setup to determine density at given altitude 

CALL SFTUP( A,ALT,T«RH0) 

0ME=To?921 1E8F4P4E~05 

VE= SORT ( ( A ( 2 ) +0ME-’'^A ( 4 ) ) •i*'*2+ ( A ( 5 ) --OME-K-A (11) -K-^a+A C 8 ) **2 ) 

FDRO=-,R-X-CDA*PHO*VF-»-1 oF+3 

COMPUTE DRAG ACCELERATION COMPONENTS 

XODno=FnPG-5(- ( A ( 2 ) +OME*A ( 4 ) ) 

YDDDPs'^DPG* (ACS) A { 1 ) ) 


7D0DP=FnRG*A ( 8 ) 
R2=R*P 


R3=PF*R 


DS=p,' 3 *pp 


P 7 =Rc^*o ;5 • 

R = A ( 7))!-A (7 ) /RF 
BSBB=9»0 

C**-«-* COMPUTE TwO-BODY ACCELERATION COMPONENTS 
X280DY=-3.9R6032F+0S-!(^A( 1 ) /R3 
Y2B0DY=-3.Q86032F+05*A( 4)/R3 
72B0DY=-3.P86032E+ns*A( 7)/R3 

C***-!!- COMPUTE EFFECTS OF EARTH, S OBLATENESS C J2 , J3 , J4 ) 
XJ20P=-2.632S1 7085572+1 0*A C 1 ) * ( 1 . -5 . ) /R5 
X J30R=+5, 94697 1 7S476E + 1 1*A ( 1 )*( 3 . -7. ) *A ( 7 ) XR7 

XJA0P=-5.1P486592772E+1S*AC 1 ) «■ ( . 42857 1 42P57-6 . -Si-B+BBeB ) /R7 
YJ20R=-2.63251 70P557E+1 0)?-A (4 )*( 1 .-S,*B)/R5 
YJ30B=+5.94697175476E+1 1 *A ( 4 ) *A ( 7) * ( 3 .-7 . *B ) /R7 
YJ40B=-5. 1 9486592772E+1 5*A(4 )* < . 42857 1 42857-6 .*B+BBBB ) 7R7 
ZJ20B=-2.63251 708557E+1 O-SS-A (7)*( 3 .-54 ) /R5 

ZJ30B=-3.5681 83053E+ 1 1* ( 1 .-1 0.*B+I 1 «6666666667*B*R ) /R5 
Z J40B=-5 . 1 9486592772E+ 1 5* A ( 7 ) * ( 2 . 1 42857 1428-10. *B+BBBB ) /R7 

C*-)?-*-)!- COMPUTE TOTAL ACCELERATION COMPONENTS 
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A { 3 ) =X3B0DV+X J30P+XJ30B+X J40B+XDDDR 
A ( 6 ) =Y^>POr>Y+Y J?OR+Y J30R+Y J40P+YDDDP 
A ( 9 ) =7;=R0DY+Z J20B+:? J30B+7 J40B+7DDDR 
903 PRTUPM 

Q0<^ NjP-MP4 ] 

IF (NHTP-NP) 904 4 904 « 903 
90 A continue 

C**** COMPUTE OSCULATING ORBITAL ELEMENTS- 
Xf<FPTH = 39RS03«? 

HX= A ( 4 I * A ( P > - A ( 7 ) * A { E ) 

HY^AC7)*A 1 I^f-ACS) 

H7=A{ 1 )*A{EI-A(A5*A(?) 

HSO= SORT ( HX-s-HX+HY-s-HY+HZ^HZ ) 

C3 = A C 7 ) * A ( 2 ) 4-A ( == ) *A ( E )+A ( 8 > *A C B ) -2 . *XKERtH/R 
C-S-^f-K-T- SEMI -major axis 
AXTS=-Xt<ERTH/C3 
XN=SORT( XKEPTH/AXIS**3) 

ECCENTRICITY 

ECCEN= 1 . +HSQ*HSQ*C3/'XKFRTH/XI<FRTH 
lE(ECCEN) 1 9«30»20 
10 ECCEN=0. 

GO TO •’0 

PO ECCEN=SQRT(FCCEN) 

I E ( ErOEN- a 0900 1 M O , 1 O , 30 
. INCLINATION 

30 FINC=ATANE CSORT{HX*Hx+HY*HY) iHZ) 
IE(EINC*nPR-.000 01 >40 »40 ,E0 
40 FlNC=Oo 

c;ETI = oTRUE. 

GO TO 70 

50 IF (E I NC^DPR- 179.99999} 80»60,60 

60 EINC=PIE 

SET I = .TRUE. 

70 CAPW=0. 

GO TO 1 on 
ascending mode 

BO CAPX\; = ATAN3(HX4 (-HY) ) 
ion TEMP=MEQ-S-HSO/XKFPTH 

true anomaly 

ANOM=R*R*ECCEN*ECCEN- ( TEMP-R ) * ( TEMP-R ) 
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IF( ANOM.GToOo ) GO TO 1;?0 
1 1 O A'N0M = 0. 

GO TO lAO 

1?0 TF(FGCFN) 101»110»101 
101 GONT.IMUF 

An‘0M = ATAN?(GQRT{ ANOM) «TEMP-R) 

TEMP= ( A ( 1 ) »A ( ? ) +A ( 4 ) * A { =; > +A < 7 ) ( 8 1 ) /R 

IF( TEf^PoLT.n, ) GO TO 130 
ANOM = ARSC ANOM ) 

GO TO 1 40 

138 ANOM=-ABS ( ANOM ) 

140 TEMP2>=A( 1 )-»COSCCAPWl+AC4)*SIN!CCAPVJ) 
TEMP=R*R-TEMPP*TEMP2 
IF( TEMP.GT.O* )GO TO 150 
TF^MD3=0. 

GO TO 190 

1 50 TFMP3= ATAN2 ( SORT ( TEMP ) * TEMP? ) 

IF( A (7) .UT.8, ) GO TO 170 
160 TEMP3= AP<; < tempo ) 

GO TO 190 

170 IF(FCCFN) 171 « ISO, 171 

171 continue 

IFCA(7)) 180,160,180 
180 teMP.3=-APS(TEMP3 ) 

190 IF(FCCEN) 250i? 00«250 
200 IF(OrTl) GO TO 240 

IF(A(T)) 200,210,209 

argument of perigee 

209 CONTINUE 
oMAlifsTFMDO 

r-Ci TO 260 

210 IF(A(8)1 220,220,230 
22-^ fMAw=o, 

GO TO 260 

23 0 SMAv<;=PIF 

GO TO 860 

24 0- 8MAw = ATAN2{A(4) ,A< 1 ) ) 

GO TO P60 

258 ?MAw=TFMP3-AN0M 
260 CONTINUE 
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wfAn anomaly 
mFAN=XN*T 

MF AN=AMOD ( MFAN » 6. ?B3 1 ) 

MFAm=MFAN*DPR 

ANOM= AMOD ( ANOM ♦ 6* aS3 1 854 ) 

SMAW= AMOD ( SMAW f 6*283 1 854 ) 

CAPw=AMOD ( CAPW t 6*283 1 854 ) 

IPO IF continue 

A MOM = A mom -K DPR 
P03 RFTURM 
FNn 

• for 4 I 5 TDAMFM , TPANFM 

subrout I NE TRANFM (AXIS, ECCEN , I NC , ASNOD , ANOM , ARGP , XS , YS , ZS , XDS , YDS , 
*R7DS) 

implicit REAL-S-8 {''-H*0-Z) 
double precision INC 

c**** THIS ROUTINE TRANSFORMS OSCULATING ORBITAL ELEMENTS TO POSITION 
C AND VELOCITY COMPONENTS IN A SPACE FIXED EPHEMER I S COORDINATE 

C SYSTEM (FARTH-CENTERED 

XKERTH =398603 * 2 
DPR=5-7, ?oc;77q<^ 

INC= INC /DPR 
A SNOD= ASNOD/DPR 
AN0M=AN0M/DPR 
APGd=ARGP/DPR 

tempi =( AXIS*( 1 *-ECCEN*ECCEN) )/( 1 *+ECCEN*COs ( ANOM ) ) 

TEMP?=tFMP1^?-S INC argp+anomj 

TFMP3=TEMP 1 *COS ( ARGP+ANOM } 

TFMP4=tEMPF*COS( INC) 

TFMP5=TEMPF*S I N t T NC ) 

TEMP6 = S I N < ASNOD ) 

TEMP7=COS< ASNOD ) 

XS=TEMP3*TEMP7-TEMP4*TEMP6 
YS=TFMP3-»TFMP6+TEMP Attempt 
7S=TFMP5 

TPMP = f 4 *XKFRTH/TFMP 1 -XKERTH/AX I s 
IFCTFMP.OT.n, )GO TO 20 
TFMPI 0=0* 

GO TO 30 

20 temp 10=SORT( TEMP) 
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30 tempi 1=ATAN?<HCCEN-«-SIN( ANOM) t 1 .+ECCEN*COS( ANOM) ) 

TEMPP= I N ( APGP+ANOM ) 

TFMPP=C0S ( APGP+AMOM ) 

TFMP?=TEMP1 0-5^ {_S‘I N( TEMPI 1 )*TFMP8+COS< TEMPI 1 ) #TEMPP ) 
TFMP3=TEMP10-«-( SIN (TEMPI 1 > *TEMP9-C0S ( TEMP 1 1 ) «"TEMP8 ) 
TFMPA=TFMP?«-C0S( INC) 

TEMdc;= tEMP2*S INC INC) 

XDS=TFMP3*TEMP7-TEMP4*TEMP6 

YDS = TEMP3-5f-TEMP6+TEMP4*TEMP7 

7DS=TFMPF 

RFTUPN 

END 

• FOR , I S BLOK , BLOK 
BLOCK data 

COMMON/D ATAB/XLABC 12) » YLABI (12) » YLAB2C 12 ) » YLAS3( 12) « YLARA C 1 2 ) » 
IYLABBC 12) »YLAB6( 12) »YLAB7( 12) i YLABS ( 12) « YLAB9 ( 12) »YLAB10{ 1 
2YLAB] 1(12) « YLABI 2(1 2) « YLABI 3 ( 12) » YLABI 4 ( 12) , YLABI 6 ( 12) »YLA 
3YLAR17( 12) »YLAB18( 12) fYLAB19( 12) *YLAB20( 12) ,YLAB21 (12) 


DATA 

XLAB 

/6H 

4 

6H 



6H 

time » 

6HIN MIN, 

6HUTES 


F 6H 


6H 

4 

6H 



6H 


4 

6H , 

6H 


F 6H 

/ 












DATA 

YLABI 

/6H 

4 

6H 

X 


6H 

IN 

KM# 

6HS » , 

6H 


S 6H 

♦ 

6H 

4 

6H 



6H 


4 

AH , 

AH 


F 6H 

/ 












data 

YLAB2 

/6H 

4 

6H 

Y 


6H 

IN 

KM# 

AHS • , 

AH 


-T 6H 


6H 

4 

6H 



6H 


# 

6H , 

AH 


F 6H 

/ 












data 

YL AB3 

/6H 

4 

6H 

7 . 


6H 

IN 

KM# 

AHS . , 

6H 


F 6H 

« 

AH 

4 

6H 



6H 


4 

6H , 

6H 


F ftH 

/ 












DATA 

YLAB4 

/6H 

4 

6H 

R 


6H 

IN 

KM* 

6HS » 1 

AH 


F 6H 

» 

6H 

4 

6H 



6H 


4 

AH , 

AH 


F AH 

/ 

• 





- 






DATA 

YLAB5 

/6H 

4 

6HX 

MIN 


6 HUS 

TW0.6HB0DY X»6H 



FAH 

♦ 6H 

46H 


<6H 



* 6H 

,6H 

,6H 


DATA 

YLAB6 

/6H 

4 

6HY 

MI 


6HNUS 

• 

6HY TWOB. 

6HODY 

4 

•F 6H 

4 

6H 

4 

6H 



6H 


4 

6H ♦ 

6H 

4 
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<p 6H •/ 

data YLAB7 /6H 
S;6H ,6H 

data YLAR8 /6H 
<F> 6H « 6H 

DATA YLAB9 /6H 
T 6H « 6H 

6H / 

DATA YLAB10/6H 
.<?; 6HFGP. « 6H 
<? 6H / 

DATA YLAB11/6H 
,6H 

S 6H Y 

data YLAB12/6H 
S 6HIN DBG« 6HS. 
<P 6H / 

DATA YLAB13/6H 

♦ 6HDFGS* < 6H 

* 6H / 

DATA YLAB14/6H 

5; 6HDFGS<- ♦ 6H 
? 6H / 

data YLAS15/6H 
S 6H ♦ 6H 

«*; 6H ■ / 

data YLAB16/6H 
6H « 6H 

<E 8H / 

DATA YLA017/6H 
T 6HGS* « 6H 
T.8H / 

DATA YLAB18/6H 
£ 6H « 6H 

S 6H / 

DATA YLAB19/6H 
« 6HDFG.S. * 6H 
£ 6H / 

DATA YLAB20/6H 


« 6H2 M3N,6HUB TW0*6HB0DY Z*6H ♦ 

«6H ,6H «6H ,6H «6H 

, 6HSTAT10. 6HN AXIS* 6H IN KM, 6HS. » 
, 6H ,6H ,6H ,6H 4 

* 6HSTATiO» 6HN ECCE , 6HNTRICI, 6HTY , 

, 6H , 6H , 6H » 6H , 

, 6HSTATIO, 6HN INCL, 6HINATI0, 6HN IN D, 
, 6H , 6H , 6H ♦ 6H « 

, 6HSTATIO, 6HN NODE, 6H IN DE, 6HGS. , 
, 6H , 6H , 6H ,6H , 

, 6HSTATIO, 6HN ARG , 6H OF PE, 6HRIGEE 

, 6H , 6H , 6H , 6H « 

9 6HSTATIO, 6HN TRuE , 6H ANOMA, 6HLY IN 9 
, 6H 9 6H 9 6H ♦ 6H 9 

, 6HSTATI0, 6HN MEAN, 6H ANOMA, 6HLY IN , 
, 6H , 6H 9 6H , 6H 9 

,6HM0DULE, 6H AXIS ,6HIN KMS , 6H. , 

, 6H , 6H , 6H 9 6H 9 

,6HMODULE, 6H ECCEN, 6HTRICIT,6HY , 

, 6H , 6H , 6H 9 SH 9 

,6HM0DULE9 6H INCLI, 6HNAT10N,6H in DE, 

9 6H , 6H , 6H 9 6H , 

,6HM0DULE, 6H node ,6HIN DEG, 6HS. , 

,, 6H ♦ 6H , 6H , 6H , 

,6HM0DuLE, 6H ARG 0,6HF PERI, 6HGEE IN, 

, 6H , 6H , 6H 9 6H , 

,6HM0DULE, 6 H true 96 HANOMAL, 6HV IN D» 
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* ♦ 6H « 6H « 6H * 6H ♦ 6H « 

♦ 6H / 

DATA YLAB 21 / 6 H ♦6HMODULE1 6H MEAN ♦ 6 HAN 0 MAL» 6HY IN D» 

‘S 6 HF^ 5 S» ♦ f)H 4 6 H . 4 6 H 4 6 H 4 6 H 4 

T6H / 

FND 

•FOR4IS TOBODY4TOBODY 

SUBROUT I NE TOBODy ( A 4 AX , 0 4 AX I S 4 S 4 X I 4 CAP , SMAw 4 TRUE 4 MEANP 4 T ) 
IMPLICIT REALMS (A-H 40 -Z) 

DIMENSION A ( 1 4 ) 4 AX{ 14 ) 4 B( 15 ) 

DOUPLF PRECISION MFANPtMFAN 
COMMON/XK/NP 4 NDTP » H 
N = 3 

C THIS ROUTINE INTEGRATES THE EQUATIONS OF MOTION USING A FOURTH 

C ORDER RUNGE-KUTTA NUMERICAL SCHEME 

1 CONTINUE 

C CALL CONIC TO GET INITIAL VALUES OF ACCELERATION 

CALL CONIC ( ] 4A4R4AXIS4ECCEN4ANOM4CAPW4SMAW4MEAN4FINC4T) 

C CALL CONIC TO COMPUTE ORBITAL ELEMENTS FOR INITIAL AND SUBSE- 

C OUENT TIMF PERIODS 

CALL CON IC<24A4R4AXIS 4ECCEN , ANOM 4 C APW 4 SMAW 4 MEAN 4 F I NC 4 T ) 

C store POSITION4VELOCITY4 ACCELERATION AND ORBITAL ELEMENTS TO BE 

C returned to main ROUTINE 

0 = R 

AX( 1 )=A( 1 > 

AX(2)=A(2) 

AX( 4 )=AC 4 ) 

AX( 5 )=A( 5 > 

AXr?) =4 (T) 

4X(F ) = /\ ( P.) 

AXIS=AXI S 

S=ECCFN 

XI=FINC 

CAP=CAPW 

SMAW=SMAW 

meanp=mean 

TRU'^^ANOM 

C COMPUTE FIRST HALF-STEP THE FIRST TIME 

C=H/ 2 e 
K = 3 *N 
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F=0-«-4, 

F=C/3. 

T=T+C 

L=1 

DO 2 1= 1 «K,3 
R(L) =A( I ) 

R (L+1 )=A( I+l ) 

B(L+?)=At 1+25 

A ( n =A( I 5+C*A ( I + n+D*A( 1+2) 

AC1 + 1)=A(I + 1 5+C*A{ 1+2 5 

2 L=L+R 

C - COMPUTF FIRST HALF-STFP SECOND TII^F 

CALL CONlC( 1 »a»r,axis«eccen»anom,capw»smav\/»mean,finc»t) 

L=1 

DO 3 1=1 ,K, 3 
8(L+3)=A( 1+2) 

A( I+l )=B(L+1 )+C*A( 1+2) 

3 L=L+5 

C COMPUTE SECOND HALF-STEP THE FIRST TIME 

CALL CONICC 1 » A«R»AXIS.ECCEN,AN0M.CAPW»SMAW»MEAN»FINC»T) 

L=1 

T=T+C 

DO 4 T= 1 «3 

BCL+4)=A( 1+2) 

A( I+l )=H-K-A{ I+2)+RCL+l ) 

A( I )=R(L)+H*B(L+1 )+F*A( I +2) 

A L=L+R 

C COMPUTE FINAL CONDITIONS FOR THE TIME STEP 

CALL CONIC ( 1 I A»R « AXIS .ECCEN» ANOM» CAPW^SiMAW^MEAN^FINC^T ) 

L=1 

DO R 1=1 *K,3 

A{ I )=B(L)+H*(R(L+1 )+F*{B(L+2)+B(L+3)+8(L+4) ) ) 

A< I+l )=B(L+1 )+F^f(B(L+2)+A( I +2 ) +2 ** ( B ( L+3 ) +D { L+4 ) ) ) 

R L=L+R 

IF(NP-NDTP) 1^37,37 
37 RETURN 

end 

•FOR, IS CONIC, CONIC 

subrout I NE CON I C ( J J , A , R » AX I S , ECCEN , ANOM , CAPw , SMAW , MEAN , F I NC , T ) 
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IMPLICIT REAL-K-8 (A-H,0-Z) 

COMMoN/XK/NP . NDTP » DT 
DOUBLE PPECISION MEAN 
LOGICAL SFTI 

THIS ROUTINE DETERMINES THE TOTAL ACCELERATION COMPONENTS AND 
ALSO COMPUTES ORB-ItAL ELEMENTS FROM POSITION AND VELOCITY 
n I MP nS I on a < 1 4 I 

SET I = , FALSE, 

PIE=T. I4IB01 
0PR = B7.E9E*77Q<=i 
I F ( J J~ 1 ) 708 , 70? , 90S 
70? R= SORT! A( 1 )**?+A(4)**?+A(7) **?) 

9S6 R2=R#R 
R3=P?-»'R 

C**** COMPUTE TWO-BODY ACCELERATION COMPONENTS 
AC3)=-S,98603?E+0E*A( D/R3 
A(6) =-3*98603?E+05*A( 4)/R3 
A{9)=-3,o8603?F4-0S-!fA( 7>/R3 ‘ 

«o? return 

908 NP=NP+1 , 

IE (NDTP-NP) 904 « 904*903 
904 CONTINUE 

C*#** COMPUTE OSCULATING ORBITAL ELEMENTS 
XKERTH = 398603,? 

HX= A ( 4 ) -«-A < B I -A ( 7 ) *A ( 8 ) 

HY=A<7)*AC2)-A( 1 )*A(8) 

HZ=A < I > *A C 8 ) -A ( 4 > -S-A ( 2 ) 

HSO=SORT( HX*HX+HY*HY+HZ*HZI 

C3 = A ( 2 ) *A { 2 W-A { 51 *A { 5 )+A ( 8 ) *A { a ) -2,*XKERTH/R 
SEMI -major AXIS 

AXIS=-XKERTH/C8 

XN=S0RT(XKEPTH/AXIS**3I 
c#*** ECCFNTPICITY 

FCCFN= 1 , +HS0*HSQ*C3/XKERTH/X!<ERTH 
IF(PCCFN) 10*30,20 

in eccfn=o# 

GO TO 30 

20 ECCEN=SORT(ECCEN) 


IF (ECCFN-. 00001 ) 10, 10,30 
C#-!;-** INCLINATION 
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30 FINC = ATAN3(30RT(HX-»-HX4-HY*HY) .HZ) 
IF(FINC*OPR-900001 )40. 40,50 
40 FINC=0, 

5FTT = .TRUF. 

GO TO TO 

50 IF( FINC*DPR-I 79.99999) 80.60.60 
60 FINC=PIF 

SFT I = .TRUF. 

70 CAPw=0* 

GO TO lOO 

C**** ASCFNDING NODF 

ao CAPW=ATAN?(HX, <-HY) > 

1 00 TFMD=HSO-«-HSQ/Xf<FPTH 
C**** TRUE ANOMALY 

ANOM=R^^■R■»•ECCEN■!^■ECCEN~ ( TEMP-R ) * ( jemp-R) 
IF{ ANOM.GT.O. ) GO TO 120 
IlO AMOM=0. 

GO TO 140 

120 IF(EGCFN) 101.110.101 
101 CONTINUE 

AN0M=ATAN2CSQRT( ANOM) .tEMP-R) 

TFMP={At 1 )*A(2)+A(4)*A(‘=^)+A{7)*A(8) ) /R 
TF( TFMD.lt. 0, ) GO TO 130 
ANOM=/^Rc (ANOM ) 

GO TO 140 

130 AM0M=-A8S<AN0M) 

140 TEMP?=A( 1 )*C0F(CAPW)+A(4)*SIN(CAPW)- 
TEMP=R*R-TEMP2*TEMP2 
IF(TEMP.GT.0. )G0 TO 150 
TFMP3=0, 

GO TO 190 . 

1 50 TEMP3=AT AN2 ( SORT { TEMP ) . TEMP2 ) 

IF( A(7) .LT.O. ) GO TO 170 
160 TFMP3= ARS ( TFMpOl 

GO TO IPO 

170 IF(ECCEM) 1 71 , 180 . 171 

171 CONTINUE 

IF(A(7)) 180. 160. IRC 

180 TFMP3=-ABS ( TEMP3 ) 

IPP IF{PCOFN) 250,200.250 
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pan IF(SPTI) GO TO ^40 

IFfA(7)) 20q,210,J?09 
•*** APGUMFMT of PFRIGFF 
POQ CONTI MUF 

c;MAv^/=T'^MP3 
GO TO 360 

31P IF(A(8)) 2?0.320«2GO 

320 SMAW=0, 

GO TO 260 
230 S<^AW=PIF 
GO TO 260 

p^O SMA«i/=aTAN? ( A ( 4 ) < A ( 1 ) ) 

GO TO 260 

2F0 G'^Aiv=TFMP3-AN0W! 

260 CONTINUF 

MFAN ANOMALY 

mfan=xn*t 

MEAm=AM0D{MEAN»6. 2831864 ) 

mean=mfan-«-opr 

ANOM=AMOD ( ANOM ♦ 6 • 283 1 864 ) 

5^MAw=AMOO f FMAW» 6.2831864 ) 

CAPw = AMOD(CAPW» 6*2831 864 ) 

■012 CONTINUF 

ANOM = ANOM * DPP 
P03 PFTUPN 
FND 

OR* IS SFTUP* SETUP 

SUBROUT I NF SFTUP < A * ALT * T I ME ♦ RHO ) 

IMPLICIT REAL*8 (A-H,0-Z) 
dimension A( 1 4 ) * UBARC 3) » FTENBI 200 > 

THIS ROUTINE SETS UP ALL ' INFORM ATI ON NECESSARY TO COMPUTE DENSITY 
IN ROUTINE JACHIA 

DATA DPR/57. 295780+0/, eCL I PT/23. 44 36D+0/*FLAG/0.0D+0/* END I D/6HEND 
s / 

IF( FLAG) 8* 10*8 
lO CONTINUE 

READ MODIFIED JULIAN DATE = JULIAN DATE - 2400000.5 
RFAD< =^* lO=:)XJO 
lOF F0RMAT(D12.P) 

read the number of INPUT DATA TO BE FOUND ON SUCCEEDING CARDS 
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RFAD('=?i ] 00)K 
ICO FORMAT < 3 T3) 
no 76 T=I »K«6 

read decimetric solar flux(ave* VALUE) .decimal year .flux .year .etc 

76 read ( 5. lOl ) FTENBC I ) .FTENBC I+l ) .FTENBI 1+2) .FTENBC 1+3) .FTENBt I +4 ) . 
*FTENPC r+F) 

101 FORMAT! 6D12.F) 

FTENP(K+1 )=FNOin 
XDAY6=XJD-36203. 

P ■ CONTI MUF 

XDAyS=XDAYS+TIME/86400. 

DAYM0=AM0D(XDAYS»365.2F) 

YFRR= 1 «5R , + ( Xn AYF/365 • 2F ) 

DETERM I NEVALUE OF SOLAR FLUX FOR DATE IN QUESTION 
M= 1 

030 IF(FTENB(M)-“ENDID) 1050. 1040. 1050 
FI0B=FTENB (M-2) 

GO TO 1 100 

050 IF(FTENB{M+1 )-YERR) 1060. 1070. lORO 
060 M=M+2 

GO TO 1030 
070 F10B=FTFNBCM) 

GO TO 1 1 00 

OBO IF(M-l ) 1070. 1070. 1090 

090 F10B=FTENB{M-2)+(YERR-FTENBCM-1 ) )*(FTENB(M)-FTENB(M-2) )/(FTENB(M+1 
*) -FTFNB tM~l ) ) 
lOO CONTINUE 

determine geomagnetic index. ap« from solar flux number 

260 IF(F10P.GE«130. ) GO TO 1280 
IF(FinB.GE.80. ) GO TO 1270 
AP=6. 

GO TO 1440 
270 aP=B« 

GO TO 1440 
280 AP=11« 

440 continue 

eio=fiob 

XLAMS= o 0 1 7203*XDA YS+. 0335*S I N( . 0 I7203*XDAYS ) - 1 .4 1 

001 IF(XLAMS-6. 2831853) 10003,10003. 10002 

002 XLAMS=XLAMS-6. 283 1853 
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GO TO 10001 
103 CONTINUE 

XL.G=CO<= (XLAM3 ) 

TFMO^SINIXLAMGI 
TFMR=C0S ( ECL I PT/OPR > 

TFMG=S I N { FCL I PT/DPR ) 

XMG=TFMP*TFWtO 

xng=tfms*tfmo 

COMPUTE right ASCENSION AND DECLINATION OF SUN 
ALSUN= ATAN2 ( XMS » XLS ) 

OLSUN=rASIN(XNS ) 

URAR( 1 1-A( 1 ) 

URAR(3)=A(4) 

URAR(3)=A(7) 

HKM=ALT 

CALL DENSITY ROUTINE 

CALL JACH I A ( HKM , F 1 0 ♦ F 1 OE , ALSUN * DLSUN < AP * UB AR . DAYN0« RHO ) 

RHO=l ,D+3*RHO 
FLAG=1 . 
return 

END 

R,IS JACHIA, JACHI A 

SUBROUT I NE UACH I a ( HKM ♦ F 1 O » F 1 OBAR , ALSUN , DLSUN * AP . UBAR , DA YNO . RHO ) 
IMPLICIT REALMS (A-H«0-2) 

COMMON FN( 5) , Y( 5) »C0E(6> 
dimension UBAR(31 «L0GRH0(4) ,SBUF(4) 

real logpho,lcrhoh 

D I mens I ON TABL ( 80 ) , TABL 1(40) 

THIS ROUTINE COMPUTES DENS I TY ( GM/CM**3 ) FROM THE MSFC MODIFIED 
JACCHIA MODEL ATMOSPHERE ( 1 967 > AS GIVEN IN — SPACE ENVIRONMENT 
CRITERIA guidelines FOR USE IN SPACE VEHICLE DEVELOPMENT ( 1969 
REVISION) .DON K« WE lOER.ED I TOR 
data TABL / 


5^-0.291 18614D+01 

•»+0«0 

D+OO.+O.O 

D+OO.+O.O 

D+OO. 

S-0 « 33835 1 060+0 1 

1+0.0 

D+OO .+0*0 

D+OO.+O.O 

D+OO. 

S-0 .405 1 0482D+0 1 

. +0 • 0 

D+OO ,+0.0 

D+OO.+O.O 

D+OO, 

S-0. 4734951 OD+0 1 

. +0 0 0 

D+OO ,+OoO 

D+OO.+O.O 

D+OO. 

£-0 o 53983993D+0 1 

.+0.0 

D+OO.+O.O 

D+OO, +0.0 

D+OO. 

S-0.59884720D+01 

1+0.0 

D+Op.+O.O 

D+OO.+O.O 

D + OO. 

$-0o65 1 43S90D+0 1 

«+o,o 

D+OO.+O.O 

D+OO.+O.O 

D+OO. 
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S-0. 70578 1550+01 , +0.0 

0+00, +0.0 

D+00,+0.0 

0+00, 

S-0. 769528200+01 ,+0.0 

D+00,+0.0 

D+00 , +0.0 

0+00, 

S 44*0,00+0/ 




data TABLI / 




S-0. 930300 130+01 ,+0,0 

0+00, +0.0 

D+00,+0.0 

0+00 , 

S-0. 100070000+02, +0.0 

0+00 .+0.0 

D+00,+0.0 

D+00, 


S-O. in609000D+0?t-0.4T6837l6D-Q6»-0*476a3T16D-06»+0#0 D+00» 

■B-Oo 1 1091000D+02 i-0.31907082D-01 «+0« 1 5483856D-0 1 » -0 . 250720980^02 « 
S-0. 1 141 10000+02,-0.246601 100-01 , +0 • 708436970-02, -0 • 597476960-03 , 
®-0. 1 16560000+02, -0.4991 53 140-02, -0.863742830-02, +0.282669070-02, 
S-0. 1 18570000+02,-0,200176240-01 ,-0.264949300-01 ,+0.660371780-02, 
S-0. 120320000+02,-0.474896430-01 , -O . 44704437O“0 1 , +0 • 104184150-01 , 
$-0. 121870000+02,-0.764255520-01 , ~0 .624866490-0 1 , +0. 139999390-01 , 
4-5<-O.OD+0/ 

EOUIVALENCF (TABL<411, 'TABLl < 1 M 

TABL ( 37 ) =-0 .84989394E+0 1 

TABL(38)=-0.0 

TABL(39) =-'^,0 

TABL ( 40 1 =-0,0 

TABL 1(37) =-0 . 1 2329000E+02 

TABLl (3B)=+n. lO565329E-0n 

TABL 1(39) =-0 .79087734E-0 1 

TABL 1(40) =+0 , 1 72758 1 OE-0 I 

COE( 1 ) =-10.48047029 

COE (2) =2.844291 123E-02 

COE ( 3 ) =-3 . 62095982 1 E~05 

COE (4 ) =2. 341 193059E-0R 

COE (5)=-7.577509214E-12 

COE (6)=9,753963073E-16 

FK=8.31432E+07 

CONL = .A2A2QA4B 

PI=3. 14 159268 

TP I =6,2831 85.36 

PI4=. 78539^2 

C0N=57, 29577951 

EP9=23. 45/CON 

F10R=E10PAR 

TOB=362.+3.6*F10P 

T0P=T0P+1 .a*(Fl n-FlOB) 

TO=TOP+( ,37+. 14-«-S IN{ TPI* (DAYNO-151 . ) /365. ) )* 
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1 (FlOR*SlNH2.-K-TPI*(DAYNO-‘=iQ.)/36‘=:») ) 

RMAG=S0RT(UBAR( 1 )*URAP( 1 ) +UBAR C ? ) *UBAR < H ) +UBAR ( 3 ) *UBAR ( 3 ) ) 
PH I = A<r, I N ( UR AR ( R ) /'-'MAG ) 

H3UN=ATAN2(URARC?) «UBAP( 1 ) ) 

10 IF(HRUN.GE.O.O) go to 20 
HRUN=HSUN+TPI 
GO TO 10 

20 IF<HRUN.LE.TPI ) GO TO 30 
H3UN=HFUN-TPI 
GO TO 20 

30 hsun=hgun-alsun 
31 IFfHPUN. GF.O, )G0 TO 32 
HSUN=HSUN+TPI 
GO TO 31 

3? IFcHPUN.LT.TPn GO TO 33 
HSUN=AMOD<HRUNiTPI ) 

33 T Au= HSUN-P 1 4+ . 209439-K-S I N ( HSUN+P 1 4 ) 

IF(TAU)270,2‘50,2r0 
250 IF(TAU-Pn2O0 i2P0»260 
260 TAU=TAU-TPI 
GO TO 250 

270 IF(TAU+PI )'280,290»290 
280 TAU=TAU+TPI 
GO TO 270 
200 CONTINUE 
29= FTA=0,'^*tPHI-nLSUN ) 

PS I = APS ( 0 . 5* ( PH I +DLSUN ) ) 

300 OOUI=2,'=: 

320 RPol=GIM(PSn 

CTAU=C0S(TAU/2.0) . 

CETA=COS<ETA) 

TF(SP«:I 1 322 « 322 « 323 

322 GPsI=0. 

GO TO 325 

323 SPRI =SP3I*^^POW 

325 IFCCTAUl 326i326» 327 

326 CTAU=o, 

GO TO '»30 

327 CTAU=CTAU**POW 

330 IF(CFTA)331 ,331 ,332 
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CfTTA = n, 

GO TO OOO 

G3? CETA=CPTA*#POW 
CONTINUE 

T=TO*( 1 . + (sSPSl+CTAU*CETA-CTAU*SPSl ) ) 

0ELT=AP+100,*( 1 .-EXP( -,08*AP ) ) 

. . , TF=T+DELT 

EPR=TF-800, 

y=EPP/ < 7«=;0 , + 1 , 7F2E-04-SERR**? 5 

FT=,02P1*FXP(-Xi^*2/2,0) 

SIGMA=GT+.D0015 
040 CONTINUF 

G=(HKM~120. ) *6476.77/ ( 6386 *77+HKM) 

TZ=TE-<TE~355. )*EXP(-8IGMA*G) 

I F ( HKN- 120.) R50 .3^2. 342 
880 TEM0D=(TE~1 l0O,0)/400c0 
S=HKM/10.0 
1=8 
S=I 

S'= (HKM/10,0)-S 
IF (8) 8801,110.8801 

2801 continue 

•IF( I .L=-. 1 > GO TO 70 
IF(I,GP.P9) go to RO 

1 = I*4_T 

50 J=I+i2 
L=1 

c compute uogrho for 4 altitudes 
no 60 K=I,J,4 

LOGRHO { L ) =TABL ( K ) +TABL ( K+ l ) *TEM0D+TABL ( K+2 ) *TEM0D*TEM0D+TABL ( K+3 ) 
l-s-TE'^On*TFiV'OD*TEMnD 
L=L+1 

60 continue 

SPUFC t 1=8* (8-1 .0)*(S-2,0) 

8RUF(2)=(8“1 .0)*(S~2.0)*(S+1 ,0) 

8RUF(3)=S*{8-2.0)*(S+1 ,0) 

SRUF ( 4 ) =S* ( S- 1 o 0 ) * ( S+ 1 . 0 ) 

C COMPUTE LOG RHO AT H KM - INTERPOLATE 

LGRHOH=( (SBUF( 1 )4{L0GRH0( 1 ) )/(-6.0) ) + 

1 < ( S8UF ( 2 ) *LOGRHO ( 2 M /2 . O ) + 
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? ( ( Snur ( T ) -K-LOGRHO c .T ))/(-?« 0 ) ) + 

.T ( ( SRUF ( 4 ) *LOGRHO ( 4 ) ) /6 . 0 ) 

C RH0=DENSITY( G/CM3) = ANTILOG OF LGRHOH 

oh0=10.**LGRH0H 
RFTUPN 

Tn IF(I.LT«0) GO TO 90 
1 = 1 
GO TO 

B'O IFU.GT.IOO) GO TO 90 
r=3GO 
GO TO FO 

90 WRITFf 6» lOniHKM 

100 FORMAT( 1H0,5H H = » E 1 5, 8 , 5X t 40HEXCEEDS ATMOSPHERE TABLE OF 0 TO 
*00 KM) 

RH0=0, 

return 

1 1 O J=I*4+1 

LGRHOH=TABL ( J ) +TABL ( J+ 1 ) *TEMOD+TABL ('J+2 > *TEMOD*TEMOD+TABL ( J+3 ) 

1 *TEMOD-«-TEMOD*TEMOD 
GO TO 6F 

34? FN(1 5 =4. OE+1 1 
■ FN(-2)=7.FE+10 

FN( 3) =7.6F+I 0 


FN{ 4) =3.4E+07 
,IF( HKIv-FOn, )3F0,34F»34F 
34F i=RP=DLOG( TF)*CONL 

FN(F) =73. 13-39. 4*^RR+5.5*ERR**2 
FN ( F ) = 1 O . C**FN ( 5 ) 

ALH=COE( 1 )+TE*(COFC2)+TE*(COEC3)+TE*(COE(4)+ 
1 TE* ( COE { F ) +TE-K-COE ( 6 ) ) ) ) ) 

N=5 

GO TO 360 
350 FN(F)=n. 

N=4 

1F(ABS(HKM-120o )- l ,0E-08) 380 ,380.360 
360 Y C 5 ) =944 ,6FFE+'^5/ ( S I GMA*FK*TF ) 

Y( 1 )=?8.*Y(F) 

Y(?)=32.*Y<=) 

Y(3)=l 6.*Y(=) 

Y(4 ) =4.*Y(F) 


10 
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)/TE 

C=( 1 .-A)/( 1 .-A-»-FXP(-Slf;MA'K'G) ) 

FN { 1 ) =FN ( 1 ) *C** ( 1 • + Y t 1 ) ) *EXP ( ~S I GMA*Y < 1 ) *G ) 

FN( p) =FN < P )*C** ( 1 . + Y(2) ) *EXP(-S IGMA*Y( 2 )*G) 

FN(P) =FN(3)*C**( 1*+YC3) )*EXPt”SIGMA*Y(3)*G) 

FN(4) =FM(4)*C**( .63+Y<4) > *EXP ( -S I GMA*Y { 4 ) ) 

IF{ ARFCHK-M-FOO. )-l .0^-08)380 *380 ,3T0 
37<^ FN<5)=FN(5)*C**( i *+ALH+YC5) ) -X-EXP ( “S I GMA*Y ( 5 ) *G ) 
"=80 5UM=0, 

no T = j ,Nj ■ 

sum=8um+fn( r ) 

390 CONTINUE 

ERR= 28 . *FN ( 1 ) +32 . *FN < 2 > + 1 6 * *FN { 3 ) +4 . *FN { 4 ) +FN {'5 ) 

RHO=ERP-«-l .66F-P4 

RFTUPM 

END 

• MAP* IM A ,R 

LIR 8YSF-XM8FCS. 

• XOT P 


764P.4=: 


• 1 

55.0 

0.0 

0.0 

0.0 


7642,655777 O 

. I 

55*0 

0.0 

0.0 

0,0 


RP 

-Op 

P0«=: 

-oi 





3 

+op 

15 

+ 02 10 





443RQ 

+05 







oi ? 








13729 

+03 

1980 

+04 13314 

+03 198025 

+04 12994 

+03 19805 

+ 04 

12690 

+ 03 

198075 

+04 12371 

+03 1981 

+04 12112 

+03 198125 

+ 04 


• FIN 
•FIN 
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